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ABSTRACT OF THE DISSERTATION 



The theory of gyrokinetic turbulence: 
A multiple-scales approach 

by 

Gabriel Galad Plunk 

Doctor of Philosophy in Physics 
University of Cahfornia, Los Angeles, 2009 
Professor Steven Cowley, Chair 



Gyrokinetics is a rich and rewarding playground to study some of the mysteries 
of modern physics - such as turbulence, universality, self-organization and dy- 
namic criticality - which are found in physical systems that are driven far from 
thermodynamic equilibrium. One such system is of particular importance, as it is 
central in the development of fusion energy ~ this system is the turbulent plasma 
found in magnetically confined fusion device. In this thesis I present work, mo- 
tivated by the quest for fusion energy, which seeks to uncover some of the inner 
workings of turbulence in magnetized plasmas. I present three projects, based on 
the work of me and my collaborators, which take a tour of different aspects and 
approaches to the gyrokinetic turbulence problem. 

I begin with the fundamental theory of gyrokinetics, and a novel formulation 
of its extension to the equations for mean-scale transport - the equations which 
must be solved to determine the performance of magnetically confined fusion 
devices. The results of this work include (1) the equations of evolution for the 
mean-scale (equilibrium) density, temperature and magnetic field of the plasma, 
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(2) a detailed Poynting's theorem for the energy balance and (3) the entropy 
balance equations. 

The second project presents gyrokinetic secondary instability theory as a 
mechanism to bring about saturation of the basic instabilities that drive gy- 
rokinetic turbulence. Emphasis is put on the ability for this analytic theory to 
predict basic properties of the nonlinear state, which can be applied to a mix- 
ing length phenomenology of transport. The results of this work include (1) an 
integral equation for the calculation of the growth rate of the fully gyrokinetic 
secondary instability with finite Larmor radius (FLR) affects included exactly, (2) 
the demonstration of the robustness of the secondary instability at fine scales {kpi 
for ion temperature gradient (ITG) turbulence and kpe <S 1 for electron temper- 
ature gradient (ETG)) which rules out the possibility that ultra-fine streamers 
could produce significant transport, (3) a demonstration that the variation in 
the phasing of the primary mode (which depend on the values of the equilib- 
rium scale lengths of the system) effects the strength of the secondary instabihty, 
distinguishing the gyrokinetic model from a previous gyrofluid model, (4) param- 
eter scans for the mean-scale gradient lengths which suggest a possible role of 
secondary instabilities in the Dimits shift and the formation of electron internal 
transport barriers (ITB) in tokamaks, (5) a formulation of the theory for fully 
gyrokinetic ions and electrons in order to explore the transition between ETG and 
ITG scales and (6) demonstrate the existence of a mechanism for the saturation 
of long-wavelength ETG modes in this ETG-ITG transition range (modes which 
have been demonstrated in simulations not to saturate when employing the ETG 
Boltzmann-ion gyrokinetic system). 

The final project is an application of the methods from inertial range under- 
standing of fluid turbulence, to describe the stationary state of fully developed 
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two-dimensional gyrokinetic turbulence. This work explores the relatively new 
idea of a phase-space cascade, whereby fine scales are nonlinearly generated in 
both position space and velocity space, and ultimately smoothed by coUisional 
entropy production. This process constitutes the thermodynamic balance which 
occurs in the true steady state of a turbulent plasma, including those found in 
fusion devices. The results of this work include (1) exact third order relations (in 
analogy to Kolmogorov's four-fifths law), (2) phenomenological scahng theories 
for the forward and inverse cascades, (3) a detailed description of the relationship 
of the two-dimensional gyrokinetic cascade to the Charney-Hasegawa-Mima and 
two-dimensional Navier-Stokes cascades, (4) a Hankel transform formalism for 
treating velocity scales in the distribution function and (4) power law predictions 
for the phase-space free energy spectra. 
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CHAPTER 1 



Introduction 

The work of this thesis addresses related but distinct problems in the theory of 
gyrokinetic (GK) turbulence, or gyrokinetics. This introduction will motivate the 
problems addressed by this thesis and outline the main results of the work. 

1.1 Gyrokinetics 

When inter-particle collisions in a plasma are too weak to maintain local thermo- 
dynamic equilibrium (over time-scales of interest) the velocity distribution of the 
particles becomes an important dynamical feature, and a kinetic description (i.e. 
the Fokker-Planck equation) is necessary to capture this. In kinetics, a plasma is 
described by a scalar field over 6 dimensions in phase space. This is clearly more 
complicated than a fluid description, over 3 spatial dimensions. The complexity 
of plasmas is furthered by its famously rich "zoology" of waves and instabilities, 
existing over an enormous dynamical and spatial range of scales. 

In magnetized plasmas, this complexity may be reduced substantially when 
the dynamics of interest occur over time-scales much longer than the cyclotron 
period, the period of rotation of particles about the mean magnetic field. This is 
the starting point for the development of gyrokinetic theory. The complexity of 
the problem is reduced in two ways. First, phase space is reduced by one dimen- 
sion (by eliminating the gyro-angle of particles about the magnetic field). Second, 
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and most important, the dynamical range of the problem is reduced by averaging 
away behavior at the cyclotron time-scale. This yields a vast improvement in the 
efficiency of numerical simulations. One estimate gives a total speedup of 10^^ 
achieved by the nonlinear gyrokinetic equation (when including the elimination 
of the plasma frequency scale, Debye and electron Larmor radius scales, and the 
ion cyclotron scale). "E^SB 

The historical milestones in the development of gyrokinetics are as follows. 
The evolution begins with the work of (ITH68P and (IRFGSp . The discovery of 
a generalized adiabatic invariant in the presence of low frequency fluctuations 
allowed the generalization of guiding center theory in the form of gryo-orbit- 
averaged kinetic equation, which became known as the linear gyrokinetic equa- 
tion. The linear theory was refined by introduction of the Catto-transformation 
and the inclusion of coUisional effects in the paper f lCT77p . The work by fITLSOp 
developed the application of gyrokinetics to the flux-surface equilibrium field 
geometry of axi-symmetric magnetically confined plasmas. Then a critical break- 
through came with the nonlinear formulation of gyrokinetics by (lFC82p . The 
hamiltonian formulation of gyrokinetics has also brought deeper understanding, 
and its history, as well as a general review of gyrokinetics, is reviewed in (IBHOTp . 



1.2 The Tokamak Transport Problem 

The motivation and context for the development of the gyrokinetic equations is in 
magnetic fusion. Although the applicability of gyrokinetics extends to astrophys- 
ical plasmas, the application on which this thesis will focus will be fusion. We 
begin this introduction with a statement of the problem of Tokamak transport. 
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1.2.1 Confinement and Fusion 



Following the presentation in (1Wes04p and (IFriOTIl . we now take a look at some 
simple considerations of confinement and fusion. The fusion reaction rate per 
unit volume is give by the product of the densities of the the two fusing species 
multiplied by an integral over the reaction cross section which accounts for the 
relative contribution from different energy levels of the bulk Maxwellian: 



= - I - ) — / a{e)eexp{ -)de 

\TT J \l / mi J mil 

= nin2 {av) (1.1) 

where fi = mim2/{mi + is the reduced mass of the fusing nuclei and a is 
the collisional cross section. To determine the conditions necessary for a fusion 
reactor, let's consider how the heating by fusion power must balance with the 
thermal losses. The fusion reaction between the tritium and deuterium nuclei 
produces an alpha particle and a neutron with 3.5 MeV and 14.1 MeV of kinetic 
energy respectively. The high-energy neutrons are not confined and pass from 
the plasma volume. The alpha particles are sufficiently confined to impart their 
energy back to the plasma. We may calculate the total power of this heating by 
multiplying the fusion rate 3? by the energy of a particles, = S.5MeV, and 
integrate over the plasma volume: 



= £a J d^Ynin2 {crv) 



S^V{av) (1.2) 
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where V is the plasma volume and we assume a constant density rii = n2 = n/2 
between lines two and three. The over-bar indicates volume-averaging. Now 
consider the overall energy balance equation in steady state. The total plasma 
energy is = 3 / d?Y riT and in steady-state, its decay due to losses is balanced 
by the externally applied heating Ph and the heating by alpha particles: 



W 

Pa + PH = — 

3n 



d'vT 



= (1.3) 

Te 

This introduces the important measure, te = {d\iaW/dt)^^l^^^, which is called 
the confinement time. The confinement time is the e-folding time of the total 
plasma energy due to thermal conduction. Ignition is achieved when the a particle 
heating matches or exceeds thermal losses. When this occurs, the external heating 
Ph is not needed to sustain the plasma energy. By setting P^ = 0, and using 



the expression 1.2 for a- heating power, we obtain an expression for a minimal 



confinement time needed for ignition: 

127" 

TITE = =- (1.4) 

To get a sense of the performance requirements of a fusion reactor, an approximate 
expression may be obtained for a constant temperature (fiat profile) between 10 
and 20 keV (from equation 1.5.5 in (lWes04p ) : 

uTte ~ 3 X lO^^m-^keV s (1.5) 
The typical densities found in present-day tokamaks vary from 10^^-10^" m~^. 
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Assuming the temperature range 10-20 keV, this gives an ignition confinement 
time of approximately r^; ~ 1.5 — 30 s. 

This timescale is very large in comparison to the micro-scales in the plasma 
such as inter-particle collisions, particle transit and bounce frequencies, linear 
and nonlinear timescales. This is the basis for the assumption of scale separation 
and suggests that a suitably defined time-average may be used as the ensemble 
average in the formal statistical treatment of fiuctuations. We will explore this 
idea in detail in chapter [2j 

The preceding simple calculation of the ignition confinement time sets a rough 
goal for the required performance of a tokamak-based reactor. There are several 
processes which enter in the determination of the confinement measure te- The 
basic mechanisms of confinement degradation, that are of chief concern, are col- 
lisional (classical) transport and turbulent (anomalous) transport. The former 
process is well-understood but unfortunately subdominant to the latter process, 
which is a subject of intense research, and a focus of this thesis. 

1.2.2 Classical transport 

The transport due to collisions in a magnetized plasma in local thermodynamic 
equilibrium is referred to as classical transport. The theory of this process is 
significantly more detailed for a toroidal plasma, so has been given the name 
neoclassical transport for that case. The classical and neoclassical theories are 
not the main subject of this thesis. However, there are some scenarios where 
anomalous transport is suppressed, and the neoclassical levels are observed (for 
example, transport barriers - see (1SGS99P ). Thus classical transport is discussed 
here briefiy as a baseline mechanism for transport and to outline some general 
considerations of scaling in experiments. 
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The diffusion of particles in tlie classical scenario occurs via random collisions. 
For this process, the scale of interest in the presence of a uniform background 
magnetic field is the Larmor radius. The random change in the particle veloc- 
ity implies a change in the the position of the center of the particle orbit, the 
gyro-center. Thus the step length for this process is the Larmor radius p and the 
time between steps is set by the inter-species coUisional frequency (like-particle 
collisions do not produce particle diffusion). This gives a classical particle dif- 
fusivity of D*^*"^ ~ ^iepf ~ i^eipl- On the other hand, the thermal diffusivity is 
caused by like-particle collisions and is given x^'~^^ ~ xT"* ~ ^iPi- The associated 
confinement time is te ~ Q-^/x! where a is the minor radius of the plasma. 

As an example of a tokamak which is approaching reactor conditions, consider 
the Joint European Torus (JET), which is capable of achieving peak density 
and temperatures Tj 28keV and n ^ 4 x lO^^m"^ (see flGor99p and also 
section 13.7 of (IFriOTP ). Using these numbers and the basic parameters of the 
machine {B = 3.6T, a = 1.25m), we may calculate the ion coUisional frequency, 
Ui = 3.5s~^, and ion Larmor radius, p = 6.7mm. The corresponding classical 
thermal diffusivity gives a confinement time 



Te ~ 2.8hrs. (1.6) 
Obviously, this greatly exceeds the minimal ignition confinement time (which. 



using equation 1.5, is te ~ 2.7s). However, it also greatly exceeds the actual 
confinement time observed, te ~ Is. 

Neoclassical transport theory can account for a small part of this discrep- 
ancy. Toroidal geometry introduces, among other features, additional timescales 
and spatial scales associated with the toroidal orbits of the particles. (We will 
not review the details here.) The net result is a substantial enhancement of 
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collisional diffusion. Witli the approximations of liigli-aspect-ratio and circu- 
lar cross-section, tlie neoclassical thermal diffusivity was derived by (1RHH72P in 
terms of the classical diffusivity to be x''^'"'' ~ O.Q8q'^{R/a)^^'^x^'"^ where (R/a) 
is the aspect ratio and q is the safety factor which characterizes the trajectory 
of the magnetic field lines as they trace out the toroidal surface. For the high- 
performance JET discharge (taking g ~ 4 and {R/a) ~ 2.4), these factors give 
an enhancement of about 40 over classical diffusion so that the confinement time 
is reduced to 

te ~ imin. (1.7) 

Clearly, neoclassical transport not nearly enough to account for the experimental 
observed confinement. The remainder of the thermal conduction is accounted for 
by the so-called anomalous transport as we shall now describe. 

1.2.3 Anomalous transport 

Micro-instabilities drive turbulent plasma flows, on scales ranging from the ion 
to the electron Larmor radius, thereby mixing plasma and enhancing transport, 
well beyond base levels set by classical diffusion. This "anomalous transport" 
has the affect of degrading heat confinement, so that fusion scientists widely 
view turbulence as an ailment of fusion devices. 

Let's examine turbulent transport in a little more detail. The anomalous con- 
tribution to transport of density and temperature can described by the equations 

^See also equation (14.126) of pViOTl) 
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/ dn 




\dt 






dT 


~dt 



-dxT 



X 



= -dxQx (1.8) 

where x is the radial direction - the direction in which the equihbrium temper- 
ature and pressure vary. The complexity and difficulty of the problem enters in 
that these fluxes must be calculated from a statistical description of the turbulent 
steady state. A standard expression for the anomalous particle flux (the formal 
theory will be worked out in chapter |2|) is: 



Tx = J d-'^j v^/i (1.9a) 
= j d\ j rf^k V;,(k)/i(-k) (1.9b) 

where the over-bar is a suitably defined averaging operator (to be specified later), 
h is the fluctuating part of the (gyrokinetic) distribution function, and v^^ is 
the turbulent drift- velocity. (Note that we have used homogeneity to obtain 



the expression 1.9b,) We may state the transport problem thus either as the 
calculation of the steady state averaged value of v^/i or the average spectrum 
v^/i. In other words, the transport problem is that of determining the stationary 
spectrum of fluctuations in fully developed plasma turbulence. 

Thus a central problem in the general theory of 'turbulence'-that is, obtaining 
a statistical description of the turbulent state which is capable of predicting 
mean-scale quantities such as drag coefficient conductivity, diffusivity, etc-is the 
essence of the anomalous transport problem. From a theoretical perspective, this 
is obviously a difficult, unsolved problem and although an exact solution has 
not been achieved, there are approaches to estimating the transport fluxes. The 
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simplest approach is the random walk mixing length estimate. This estimate 
assumes that the transport spectrum is dominated by a single characteristic scale 
and assumes a diffusive random-walk process such that the ffuxes are proportional 
to a diffusivity D: 



Qx ~ -nx d^T (1.10a) 

and 

(1.11) 

where Ax is the characteristic step length and At is the typical time between 
steps. If one assumes that the step length Ax is the ion Larmor radius (i.e. that 
the transport spectrum is dominated by A; ~ p~^, as is supported by nonlin- 
ear simulations) and the step time is the inverse of the linear diamagnetic drift 
frequency (c<j^ ~ (kp) Vth/Lx ~ Vth/Lx), one obtains gyro-Bohm diffusivity 



Dgb = pfvth/LT (1.12) 

where L^^ = —d\nT/dx. Assuming the high performance JET parameters from 
above (so Vth ~ 1-2 x lO^m/s and ~ Im), we obtain a resulting confinement 
time te ~ a? /Dgb ~ 3s. Clearly, this is much closer than the classical or 



neoclassical estimates (1.6 and 1.7) to the observed confinement time te ~ Is. 

The gyro-Bohm estimate for the turbulent diffusivity reveals the scaling of the 
transport fluxes (under the assumption of scale-separation) but does not give a 
specific quantitative prediction or include details of the nonlinear dynamics. How- 
ever, the validity of gyro-Bohm scaling has been tested via 

simulations«M23EWDol 
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and it is the accepted scaling in the hmit of strong scale separation between the 
Larmor radius scale and the equilibrium variation scale - i.e. the limit of small 
= p/Lt. 

As a general note on scaling, let's compare the neoclassical confinement time 
to the gyro-Bohm confinement time. In particular, the temperature dependence 
of the collisional frequency goes as T~^/^ so that the classical and neoclassical con- 
finement times will increase as T^/^ and increase with the system size as roughly 

(assuming a constant aspect ratio). On the other hand the gyro-Bohm diffu- 
sivity increases as T^/^ due to increased Larmor radius and linear diamagnetic 
frequency. This means gyro-Bohm confinement worsens at higher temperature 
(of course, the positive scaling of confinement in system size is the dominant effect 
in scaling to reactor designs). Thus the hotter temperatures expected in reactor 
conditions should make classical and neoclassical contributions to transport even 
less important, compared with the anomalous contribution. 

1.2.4 Weak turbulence 

Quasi-linear theory gives an estimate for the anomalous transport fluxes under 
the assumption of weak turbulence. 'E^^'^^'^^® This estimate is more sophis- 
ticated than the dimensional analysis of random walk estimate (given above) in 
that it includes linear mode solutions for a specific set of plasma parameters and 
equilibrium geometry. Thus it is capable of capturing the parametric dependence 
of specific features such as the anomalous particle pinch. ''^ElQll n ^jgo has the 
advantage of being built upon rigorous mathematical footing via a perturbative 
expansion in the fluctuation amplitude. However, it is limited since a linear mode 
spectrum may not adequately describe the fully developed state. The quasilinear 
theory must also provide a guess for amplitudes of the modes - the "saturation 
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amplitude." To sketch an example, let's return to the expression |1.9b and ap- 
proximate h quasilinearly. From the gyrokinetic equation (Eqn. 2.55) we may 
write (neglecting the nonlinearity and the equilibrium curvature and V-B drifts) 



- icj(k)/i(k) ~ ■ v^(k)Fo (1.13) 

where k^ = V In Fq ~ xL-^^ - assuming that the density gradient is weaker than 
the temperature gradient, i.e. L^^ ^ L~^. Also, . (Note, we have only included 
the "irreversible" contribution to h in going from the first to second line.''^^^) 



Substituting this into the expression 1.9b we get 



r. d\ j d3kx-v^(k)-^k,-v^(-k)Fo 



I' 



d^k^|x-v,(k)|2^ (1.14) 



where 7 = Im[c(j] and the second line is obtained by noting that the is real so we 
may take the real part of the expression. We may eliminate in favor of the fluid 
displacement variable, S,-, which is the time integral of the drift velocity: d$,/dt = 
v^. Thus we may write |x ■ Vp(.(k)p ~ |u;p|^^(k)p leading to the expression 



no 



r.~y d'k^j^m^ (1.15) 

Now if we assume that the dominant part of the spectrum comes from k ^ 
and furthermore assume that ^ ~ p, we obtain the quasi-linear estimate for the 
transport fluxj^ 



^The additional approximation 7 ~ vth/Lr applied to equation 1.16 gives gyro-Bohm dif- 
fusivity argued above by random walk. Thus the random walk mixing length estimate is 
sometimes referred to as a quasi- linear estimate. 
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(1.16) 

One limitation of the assumption ^ ~ fcj^ is that the typical fluid displacement 
may in fact be less than /cj^ if the saturation amplitude is low. Thus ^ ~ p may 
be considered an upper bound imposed by the geometry of the flow for the case 
where the spectrum is dominated by /c^; ~ p^^. This case is represented in flgure 



1.1(a) If, on the other hand, there are anisotropic radially elongated streamers, 
as can be the case with electron temperature gradient (ETG) driven turbulence, 
the dominant wavenumber may correspond to ky but k^. <^ ky so that we 



may have ^ P ^ see flgure 1.1(b) It is clear that the proper formulation of 
mixing length phenomenology requires insight into the fully developed turbulent 
state - a state which may or may not be successfully treated with quasi-linear 
theory. 

1.2.5 Strong turbulence 

In the strong-turbulence limit, the nonlinearity is treated non-perturbatively and 
the spectrum of turbulent excitations cannot a priori be assumed to retain lin- 
ear mode features. This problem is treated most commonly by direct numerical 
simulation. Although there is great progress in the ability to simulate turbulent 
transport, it remains a problem that pushes the limits of computational tech- 
nology and progress in theoretical understanding remains critically important. 
Strong turbulent transport may be divided into two categories: Inertial range 
turbulence and large scale or energy-containing range turbulence (see for 
instance flCG85D l 

The inertial range understanding of turbulence, pioneered by Kolmogorov'^ ^"^^^ -^ | Ko141c | [ Koi4ia |i 
is characterized by the nonlinear cascade, whereby dynamically invariant quanti- 
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(a) Isotropic turbulent mixing: gyro-Bohm diffusivity assumes ^ ky ^ 
p"^ and ^ -K/kx 




(b) Mixing by radially elongated streamer: Enhancement of transport 
by anisotropic turbulent mixing. 

Figure 1.1: Transport by turbulent mixing: The radial mixing length is 
set by the strength of the flow (the saturation amplitude) but limited by the 
dominant wavenumber in the transport spectrum < irk"^. 
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ties are transferee! locally in fc-space, to different scales, producing a self-similar 
spectrum which is independent of the specific way that the system is driven (uni- 
versality). It is a common point of view that the inertial range understanding 
of fluid turbulence has a limited applicability to the tokamak transport problem. 
Indeed, simulations suggest that the majority of transport occurs at linearly un- 
stable scales where the majority of free energy injection occurs - such scales, 
by definition, are outside of the inertial rangej^ Also there is evidence of self- 
organization in plasma dynamics at scales that are not well-separated from the 
system scale''^^ so that the traditional assumptions of homogeneity will not ap- 
ply to these cases. For these reasons, the application of inertial range concepts 
to gyrokinetic turbulence has been largely unexplored, although much is known 
about fluid models such as the Charney-Hasegawa-Mima (CHM) turbulence. 
The final chapter, chapter |4| of this thesis is devoted to an inertial range theory 
of gyrokinetics and we will argue that this theory has a place in the study of 
anomalous transport. For instance, the inertial range understanding is impor- 
tant in determining the proper resolution needed in numerical simulations. That 
is, the fine-scale inertial range transfer of energy to the coUisional scale is the 
process by which the true steady state is achieved in a turbulent plasma; and the 
coUisional scale determines the smallest scale which must be resolved in a simula- 
tion. This will be discussed more in chapter |4} At this point, we emphasize that 
although an inertial range theory can not come close to fully describing tokamak 
turbulence, we believe it is a key part of the picture and that, in general, its 
applicability is an open question. 

The range where energy injection is non-negligible is the "energy-containing" 
range because the fully developed steady-state tends to have the majority of en- 

■^It is also possible to sustain turbulence without unstable linear modes via a nonlinear 
process which draws from the background gradient in free energy. Thus, linear instability 
is not a definitive criteria that separates inertial and non-inertial ranges. 
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ergy contained here. A statistical description of the turbulence spectrum in this 
range depends on the particular features of instability drives and other features of 
equilibrium plasma which can strongly affect the nonlinear physics. It therefore 
does not exhibit the universality that characterizes the inertial range. For this 
reason, there does not exist as much of a unifying framework to describe phe- 
nomena of the energy-containing range. Theoretical progress in this subject is 
driven by (1) the intuition built by experiment and direct numerical simulations, 
and (2) theoretical models of nonlinear processes inspired, at least in part, by 
these observations. The material in chapter [3} which is mostly concerned with 
the energy containing range, falls into the second category. 

As a final note about strong turbulence, note that although the random walk 



mixing length diffusivity estimate 1.12 may be obtained within quasilinear frame- 
work, mixing length phenomenology in general does not assume weak-turbulence 
so it can be applied to strong turbulence regimes as well. This will be discussed 
in chapter |3| 



1.3 Results 

We now turn to the results of this thesis. There are three chapters in the body of 
the thesis, corresponding to the three separate but related projects undertaken 
by the author and his collaborators. These projects all investigate theoretical 
aspects of gyrokinetic turbulence and are closely related but ultimately, stand 
alone individually. The first project, presented in chapter [2] and the appendix |A} is 
a yet unpublished study of the equations of gyrokinetics extended to describe the 
mean-scale evolution of a turbulent plasma. The final two chapters are adapted, 
with little modification, from papers by the author and his collaborators - the first 
published in Physics of Plasmas a^fj the second submitted to the Journal 



15 



of Fluid Mechanics^^^ . These two works are a closer look at the details of 
fully developed gyrokinetic turbulence and are described below. 

1.3.1 Chapter [2j Gyrokinetics as a transport theory 

The theory of transport in magnetically confined fusion plasmas began with the 
formulation of classical transport theory ''^^^^ and then neoclassical transport 
theory (a review is given by two of the pioneers in (1HH76P ). These theories 
followed the example set by the theory of collisional transport in neutral gases 
pioneered by Chapman and Enskog. While being great achievements, neither 
of the classical transport theories - which considered transport by diffusive colli- 
sional processes - could, in general, come close to predicting the level of transport 
in tokamaks. 

As we have been discussing, it has become the popular belief that the equation 
of magnetized turbulence, namely the gyrokinetic equation, must be solved to 
predict "anomalous" transport levels. On the other hand, it is observed that 
under some conditions, neoclassical transport does correctly predict the base 
level of transport - i.e. ion thermal transport in transport barrier regions. Thus, 
the most complete descriptions unify neoclassical and anomalous transport. We 
provide such a description in this chapter, chapter |2} 

1.3.1.1 Scale separation, locality and gyro-Bohm scaling 

Our treatment of gyrokinetics is standard in many regards. In particular, the or- 
dering of terms and the final gyrokinetic equation obtained is consistent with stan- 
dard non-linear delta-f gyrokinetics flFC82p . The main feature which distinguishes 
our approach from others is (1) the assumption of separation of scales between 
equilibrium scales and the fluctuation scales and (2) the method (intermediate- 
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scale spatial averaging and time averaging) by which we separate fine scale fiuc- 
tuations and mean-scale quantities. This permits the simultaneous derivation of 
(neo) classical and anomalous transport theory. 

We can view separation of scale as an assumption of local homogeneity. Con- 
sider a sub-system which extends over a volume much smaller than the equilib- 
rium scale but much larger than the turbulent scale. Although there are turbulent 
fiuctuations locally, we assume that these fiuctuations do not accumulate coher- 
ently to produce variation on the size of the sub-system-for instance, the root- 
mean-square of the fiuctuation can only change significantly on the mean scale, so 
is approximately constant (homogeneous) in the intermediate-scale sub-system. 

At the level of direct numerical simulation, scale separation implies that small 
sub-domains of the full system may be simulated separately. Such a domain may 
be chosen to be a full annulus segment comprised of the full volume contained 



between two toroidal surfaces (see section A.2.3), or even smaller domains called 
fiux tubes, which cover a small volume following a magnetic field line. These 



domains are represented in figure 1.2 A reduced domain yields a large improve- 
ment in computational efficiency and has made possible the direct simulation 
of the macro-evolution of the plasma. This has been achieved in the work by 
(IBarOSP which, has used the formulation of transport (based upon the separation 
of scales) which will be presented in this chapter. 



1.3.1.2 Chapter [2] overview 

The presentation given in section[2]gives a detailed description of the assumptions 
and methods leading to gyrokinetic theory and the mean-scale transport theory - 
the solutions to which constitutes the ultimate goal of the study of turbulence in 
fusion plasmas. This material will constitute the fundamental framework upon 
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Figure 1.2: Flux tube simulation domain: A thin flux tube volume surrounding 
an equilibrium fleld line is rendered in green. Also, a small segment of the full 
annulus volume is included on the left side of the torus. A larger segment of an 



annulus volume is pictured in flgure A.l (courtesy of G. Stantchev) 
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which the subsequent chapters are built. The notation will largely carry through, 
with the exception of specific normalizations which will be introduced explicitly 
before they are used. 

The chapter will begin with introductory material describing the ordering 
assumptions, notations and then proceed to the derivation of the gyrokinetic 
equation. The ordered expansion of the Fokker-Planck equation is extended to 
one order beyond the gyrokinetic equation to obtain the transport equations for 
density and temperature. Finally a detailed look at Poynting's theorem and 
entropy balance will give a perspective of the overall energy balance and mean- 
scale thermodynamics. 

An extension of this work to the case of toroidal magnetic field geometry is 
given in the appendix |Aj The material in chapter |2] and appendix |A] constitute 
a body of work that is still developing. In light of its special applicability to the 
ITER design, it may be adapted to an ITER-specific theory. 'E^^^ However, in 
its present form it does have several complete results, some of which have already 
been incorporated into the work by (IBarOSp as mentioned above. 

1.3.2 Chapter [3f Primary and secondary mode theory 

To understand micro-turbulent transport, one typically starts with the linearly 
unstable modes which drive the turbulence. These "primary" modes in general 
exploit a gradient in the free energy of the background plasma and magnetic field. 
As these modes grow, they begin to interact nonlinearly, causing the transition to 
a turbulent state - a state where transport by the mixing of plasma would bring 
about the relaxation of the driving gradient, if an external drive were absent. In 
the case of tokamak turbulence, however, the steady state operation corresponds 
to a balance between the overall loss by transport and injection of heat and par- 
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tides. This fully developed state remains highly turbulent throughout operation 
of the device. 

This chapter will discuss some details of primary mode analysis, but is mostly 
focused on "secondary instability" theory. A secondary instability is an exponen- 
tially growing mode induced by the presence of a linearly unstable mode which 
has itself grown to sufficient-amplitude. Secondary modes or 'secondaries' are a 
key part of the nonlinear physics. Their study in the context of plasma turbu- 
lence can bring understanding of the transition to turbulence. Also, and perhaps 
surprisingly, secondary instability theory may also be used to determine features 
of the fully developed state. 

In fully developed plasma turbulence, it cannot be assumed that the chaotic 
state will bear resemblance to the unstable linear modes from which the state 
originated. However, it has been observed that radial streamers (structures ex- 
isting in the fully developed turbulent state) bear a close resemblance to unstable 
linear modes. lEsiosj LJenogli hght of this observation, it is not surprising that the 
theory of secondary instabilities has been used successfully to predict saturation 
amplitudes in fully developed turbulence. ''^^Qli These saturation amplitudes, as 
argued above, place an upper bound on radial fluid displacement and therefore 
enter directly in the mixing length phenomenology of transport. 

1.3.2.1 Chapter |3] overview 

Chapter [3] is largely taken without modification from the work (lPlu07p . As men- 
tioned, we have added some more detail on the solution of the linear dispersion 
relation. The main subject of this chapter, however, is a fully gyrokinetic sec- 
ondary instability theory. 

Our investigation into the fully gyrokinetic secondary was first motivated, in 
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part, by the possibility of an ultra-thin streamer, with kyp^ significantly greater 
than 1. Using the mixing length ideas, such a streamer could enhance transport 
via its very rapid turnover and large step length. Also, these structures could 
evade detection, as (1) the most advanced experimental techniques for measuring 
electron temperature/density fluctuations''^^^® are able to resolve k ~ lcm~^ 
(which is much larger than the electron Larmor radius for tokamak plasmas) 
and (2) the computational requirements of a simulation would be quite large to 
resolve the temporal and spatial scale range of the nonlinear state in which these 
structures exist. 

The first result coming from this work is that the very fine-scale {kpe > 1) 
ETG secondary instability is robustly unstable, over the full range of parameters 
investigated. This indicates that ultra-thin streamers are unlikely to exist. Then, 
some unexpected results are reported. By scanning across several parameters, we 
find features which suggest the role of secondary instabilities in affecting such 
macroscopic phenomena as profile stiffness and the electron transport barrier. 
Specifically, it is found that the secondary instability for the two dimensional 
ETG mode ( "toroidal" mode) exhibits a significant sensitivity to the the critical 
gradient parameter {R/LT)cv\t (which determines the instability of the primary 
mode). This finding is used, again with the aid of mixing length phenomenology, 
to offer an explaination of observed transport reduction near marginal stability 
of the primary mode (the essential idea behind the "Dimits shift") and the fact 
that ETG turbulence exhibits less profile st iff nes^ than ITG turbulence. 

Using gyrokinetic secondary instability theory, we also are able to explore the 

^From experiments and simulations of tokamaks, there is a large amount of evidence that 
the equilibrium profile is maintained at a near-marginal state by the underlying transport. This 
concept is known as profile stiffness because transport fiuxes increase sharply if the equilibrium 
is pushed beyond a critical level - which means the profile exhibits a stiffness or resistance 
against being externally tuned. 
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intermediate scale range between the electron Larmor radius and the ion Larmor 
radius. In this range, which is an example of what we later call the nonlinear 
phase-mixing range, the ion gyrokinetic equation is strongly affected by the gyro- 
averaging which cannot be approximated by a fluid treatment. We calculate the 
primary and secondary instability using a full gyrokinetic description for both ions 
and electrons and are able to locate the wavenumber of minimal secondary growth 
rate. It is argued that the trough in the secondary growth rate curve may, under 
suitable conditions, correspond to the peak of the turbulent spectrum, which is 
one of the most important features in applying mixing length phenomenology to 
predicting the diffusivity. 

To complete the study, we also investigate the secondary instability theory 
for a three dimensional slab configuration. We confirm the findings of previous 
works which indicated that the nonlinear behavior of ITG and ETG do not differ 
as strongly as they do for the toroidal case. 

1.3.3 Chapter [4j The phase space cascade in two dimensional gyroki- 
netics 

While chapter [3] is concerned with the energy containing range, chapter |4] will 
shift to a study of inertial range gyrokinetics. This is a very different approach 
to describing fully developed turbulent state, and although the inertial range 
represents a relatively small contribution to the "bottom-line" of fusion theory - 
e.g. turbulent transport - it is surely host to important nonlinear phenomena at 
play in a general turbulent plasma state. 

Gyrokinetics assumes anisotropy in the fluctuations such that structures are 
elongated in the direction of the equilibrium magnetic field, i.e. that /cy <^ k±. 
In a subsidiary limit, when the effect of fcy is subdominant to the nonlinearity. 
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we may neglect fcy altogether to obtain two-dimensional gyrokinetics. This two 
dimensional assumption is the basis of the successful Charney-Hasegawa-Mima 
(CHM) fluid model'EEiZlllHMZi which is closely related to the equation of in- 
compressible two-dimensional fluid turbulence. There are at least two scenarios, 
relevant to tokamak turbulence, in which two dimensional gyrokinetics may be 
an appropriate description. 

First consider the case of the cascade of ion-scale turbulence to scales much 
finer than the ion Larmor radius, kpi ^ l|^ First, for kp ^ 1 the dynamics 
are dominated by a nonlinear phase-mixing process which precludes a fluid ap- 
proximation. Also, in this limit the linear modes may be stable (a necessary 
but not sufficient condition for inertial range treatment) and the strength of the 
nonlinearity grows with k± and will dominate over the parallel compressibility 
term if we assume fcy is bounded. For instance, one conventional assumption is 
^li ~ 1/qR. That is, the parallel wavelength is assumed to be roughly the dis- 
tance between the inboard and outboard sides of the toroidal magnetic surface, 
i.e. the distance between the good and bad curvature regions. Another hypoth- 



esis, given by (1SCD08I1 . is that the parallel wavenumber is determined by the 
distance that particles stream along the field lines during a nonlinear turnover 
time - this implies that the parallel term remains balanced with the nonlinearity 
and cannot be formally neglected. However, it is pointed out in (ISCDOSp that 
this term becomes much less efficient as a phase-mixing mechanism and does not 
break the conservation of the forward cascading invariant. Thus this hypothesis 
of balance is not incompatible with the the exact two-dimensional theory. 

Another scenario where the high-A; two dimensional theory of gyrokinetics 
may be useful is in describing the inverse cascade from turbulence driven at the 
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Alternatively, we may consider the cascade from the electron-scale turbulence to even finer 



scales. This case relates closely to the ion-scale cascade, as is described in detail in section 4.2 
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electron Larmor radius to scales as large as the ion Larmor radiusjj (This may 
be relevant in cases where the ion scale turbulent drive is quenched by shear flow 
in, for example, a transport barrier - although further investigation is necessary 
to establish this.) This scenario is also inextricably kinetic, as the nonlinear 
phase-mixing in the ion gyrokinetic equation is strong at kpi > 1. 

1.3.3.1 The inertial range cascade in tokamak turbulence 

Notably, there are parameter regimes for tokamaks which exhibit linear insta- 
bility across nearly the full range of dynamical interest. For instance, the ITG 
instability may blend continuously into the trapped electron mode (TEM) which 
is continued as the ETG instability. Here, it is not clear that the inertial theory 
is directly applicable. However, the rigorous justification for inertial range treat- 
ment requires a demonstration that the non-conservative terms (i.e. the linear 
drive terms) become subdominant to the other terms in the dynamical equation. 
This is possible even when a linear instability is present - the linear instability 
would simply be ineffective at injecting energy. 

Even when the case of interest does not exhibit an inertial range, the classical 
inertial range cascade should be understood as one of the possible underlying 
nonlinear processes. Also, because an inertial range cascade is a universal fea- 
ture which is insensitive to the specific forcing mechanism, the verification of 
theoretical features of the cascade may be used in diagnosing and benchmarking 
numerical simulations and determine resolution requirements. 

As a final note of motivation, the fine-scale contribution to transport fluxes 

®As with neutral fluid turbulence (and CHM turbulence), the two dimensional scenario in 
gyrokinetic turbulence two dynamically invariant quantities. The spectral scaling of the two 
invariants has a fixed relationship, as with the case of fluid turbulence, which forces a dual 
cascade. 
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may be obtained from an inertial range theory. Although the majority of trans- 
port is beheved to originate from the hnearly unstable scales kpi < 1, the turbu- 
lent transport at fine scales may represent a non- negligible portion of transport. 
A forward-thinking theory of plasma turbulence should include such contributions 
as they will ultimately be part of quantitative predictions. 

1.3.3.2 Chapter |4] overview 

These preceding considerations suggest the need for a high-A; inertial range theory 
of gyrokinetics in the two-dimensional limit. We present such a theory in the 
final work of this thesis. We offer that the case of two dimensional gyrokinetic 
turbulence is a simple paradigm for kinetic plasma turbulence. 

We study the inertial range dual cascade, assuming a localised random forcing. 
This cascade occurs in phase-space (two dimensions in position-space plus one 
dimension in velocity-space) via the nonlinear phase-mixing process, at scales 
smaller than the Larmor radius. In this "nonlinear phase-mixing range," we show 
that the turbulence is self-similar and exhibits power law spectra in position and 
velocity-space. The velocity-space spectrum is treated via a Hankel-transform 
which fits naturally into the mathematical framework of gyrokinetics. We derive 
the exact relations for third order structure functions, in analogy to Kolmogorov's 
four-fifths law. 

The two dimensional gyrokinetic system bears some resemblance to the equa- 
tions of incompressible fluid turbulence and, notably, may be rigorously reduced 
(in the appropriate long wavelength limit) to the familiar Hasegawa-Mima equa- 
tion or the vorticity equation for two dimensional Euler turbulence. We inves- 
tigate the relationship between these theories. First we review the derivation 
of the CHM/vorticity equation from gyrokinetics. Then we derive a relation- 
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ship between the fluid and gyrokinetic invariants and, after a phenomenological 
derivation of the inertial range scahng laws, present a picture of the full range of 
cascades from the fluid range to the fully kinetic range. 
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CHAPTER 2 
Gyrokinetics as a transport theory 



2.1 Introduction 

This chapter will introduce the basic equations of gyrokinetics and gyrokinetic 
transport. On a basic level, the derivation is a rigorous asymptotic expansion of 
the Fokker-Planck equation centered around the assumption of scale separation 
between equilibrium and fluctuation quantities. The ultimate goals are (1) to 
obtain the gyrokinetic equation to describe turbulence at the ion Larmor radius 
scale, (2) obtain the classical counterpart to describe perturbations that vary 
smoothly in space and time and (3) obtain the evolution equations for the mean- 
scale quantities (temperature, density, magnetic field) which are influenced by 
both turbulent fluctuations and coUisional effects associated with the inhomo- 
geneity of the equilibrium (classical and neoclassical transport theory). 

For the sake of cohesiveness with the other chapters in this thesis, we will avoid 
discussion of axi-symmetric geometry, flux coordinates and neoclassical theory in 
this chapter. Although the full-geometry toroidal case has been treated (and 
is included in the Appendix), the remainder of the chapters in this thesis will 
take a local slab (constant curvature) approximation and we find it to be more 
simple and illuminating in this chapter to formulate gyrokinetics and transport 
in a local, geometry-free manner. 
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2.2 Scale separation 

The key feature which distinguishes the following approach to gyrokinetics and 
transport is the assumption of multiple disparate scales in both time and space. 
Although there are several approaches to deriving the gyorkinetic system of 
equations and multiple forms of the final equations (e.g. conservative, non- 
conservative, gyro-center density, real space density) they share a common set 
of ordering assumptions. In the following approach we take scale separation as 
the fundamental assumption and the ordering assumptions of spatial and time 
derivatives follow from this. 

From this approach, classical (or neoclassical) theory is recovered from mean- 
scale behavior (governed by the kinetic equations under a suitably defined smooth- 
ing operation) while turbulent theory manifests itself in the fine-scale behavior of 
the perturbative fields. In other words, by starting with the assumption that the 
physical fields have dependence concurrently on multiple scales, we can recover 
both theories on the same footing. Other works having achieved this synthesis of 
anomalous and classical theories are ( IShaSSt IBal90( ISH95t ISOH96all . 

It should be pointed out that this approach does have some weakness. In 
particular it is not designed to treat meso-scale phenomena''^^ where variation 
on a scale intermediate to the micro- and macro-scales is present. This is pos- 
sible with transport barriers (where the equilibrium background varies on scales 
approaching the turbulence scale). However we reiterate that scale separation 
is widely observed in tokamak experiments and should become an excellent ap- 
proximation for ITER where p* will be quite small. Thus while there are surely 
many interesting and undiscovered meso-scale phenomena, it is the view of the 
authors that the case of scale separation is crucially important to understanding 
tokamak turbulence and transport. 
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In this section we will introduce the two characteristic spatial scales (whose 
ratio defines the fundamental small parameter of gyro-kinetic theory) as well as 
the three characteristic time scales. We then discuss the method of multiple 
scales as a foundation to the derivations of this paper. Finally we detail the 
gyrokinetic ordering scheme. 

2.2.1 Characteristic scales 

We assume the existence of two important spatial scales, the macroscopic spatial 
scale length L and the microscopic spatial scale length pi. The macroscopic scale 
length is the distance over which equilibrium quantities (such as density and 
pressure) vary and is taken to be the size of the plasma-the minor radius, the 
major radius, etc., will all be taken to be order L. The microscopic scale length is 
taken to be the ion Larmor radius (we will omit the subscript unless there is need 
to distinguish it from the electron Larmor radius). These scale lengths define the 
fundamental expansion parameter: 

e = I « 1 (2.1) 

(Which is what is commonly referred to as p*.) For example in ITER these 
lengths will be approximately: L ~ n/\S/n\ ~ 2m (the minor radius of ITER) 
and p ~ 2mm^^. This gives a strong expansion parameter: e ~ 10 ^. There 
are three basic time scales of interest defined by characteristic frequencies. The 
first is the fast ion cyclotron frequency fioi (we will omit the subscript here also). 

1^0 = ^ (2.2) 

rrii 

In ITER fio ~ 2 X l{firad/ s. Next is the medium turbulent frequency uo. 
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where Vthi is the ion thermal speed. On ITER u ~ xlO^rad/ s. The third time 
scale is the slow transport time r. On ITER r ~ 3.7s. These frequency scales 
have a simple relationship in terms of e. 



r 



-1 



(2.4) 



2.2.2 Method of multiple scales 

The method of multiple scales (see for instance (1B078P ) is characterized by the 
formal introduction of independent variables to describe dependence on different 
scales. The equation describing the system of interest is then perturbatively 
expanded and solved in these independent variables. 

For our purposes, the distribution function, /(r,t), has dependence on space 
and time (and velocity). Here, however, we formally assume 



with tg = t { "slow time" ) and rg = e r ( "slow space" ) taken to be independent 
variables. (The fast time scale dependence will only enter at the level of single 
particle dyanics. In gyrokinetics, this behavior is not included in the collective 
dynamics captured by the distribution function and associated fields.) Now we 
expand the distribution function, electric field and magnetic field in powers of e 



/ = /(rs,r,ts,t) 



(2.5) 



f = Fo + 6h + 6h 



(2.6) 
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and 



B = Bo + (5B, E = 6E (2.7) 

with(5/„/Fo ~ C(e"), |5B|/|Bo| ~ C(e), etc. Additionally we take \5E\/\vthiBo\ ~ 
(9(e) where Vthi is the ion thermal speed which will be taken to be \fTifrni. This 
expansion is the starting point for the method of multiple scales. Now we must 
specify the assumptions of scale dependence for the problem of interest, gyro- 
kinetic transport theory. We assume the equilibrium quantities Fq and Bq have 
dependence only on the slow (transport) time variable tg, the slow (macroscopic) 
spatial vector rg and velocity v. The perturbative quantities 5B, 5E and bj have 
additional dependence on the turbulent scale variables denoted simply by t and 
r. Furthermore, as with standard gyro-kinetic ordering, the microscopic spatial 
dependence is assumed only in the direction perpendicular to the equilibrium 
magnetic field Bq. We can summarize these statements as follows: 



Fo = Fo(rs, v,t,) 
Bo = Bo(rs, v,ts) 
= 5/(r±,rs, v,t,tg) 
SE = 5E(rx,rs,v,t,t,) 
5B = 6B{r±,rs,v,t,ts) 



(2.8) 



These assumptions are important for all that follows. However, we find it pos- 
sible to avoid the cumbersome subscript notation entirely. This is done by first 
explicitly stating the ordering of spatial and temporal variations (described in 
section 2.3) and by introducing spatial and temporal averages that serve to sepa- 



rate scale dependence (section 2.4). In the following analysis, we thus need only 
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refer to a single space vector r and time variable t. As a final note, velocity space 
dependence is here assumed to occur on the scale of the ion thermal velocity. 
There are subtleties involved in the formation of velocity structure and these will 
be explored in Chapter |4] of this thesis. 

2.3 Gyrokinetic ordering 

Here we describe the ordering of spatial and temporal variation that results di- 
rectly from our assumptions of scale dependence. In brief, the slow scale de- 
pendence of equilibrium quantities results in space and time derivatives that are 
an order smaller than the same derivatives performed on perturbative quanti- 
ties. These orderings are essential for solving the Fokker-Planck equation and 
are summarized in the following table. 

Ordering of Spatial and Temporal Variation 

• V ~ acting on Fq or Bq 

• ^ ~ 0{e^Lo) acting on Fq or Bq 

• ~ acting on 5B or (5E 

• V|| ~ 0{{) acting on 5fi, 5B or 5E 

• ^ ~ C(co') acting on (5B or 5E 

where V_l and Vii denote the gradient perpendicular and parallel to Bq. 
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2.3.1 Electron and ion orderings 



Because p is taken to be tlie ion Larmor radius, we are considering ion scale 
turbulence. We may imagine, for instance, that the electron-larmor-radius-scale 
instability (ETG) is not present in the system of interest. (In practice, it is be- 
lieved that the feedback of electron scale turbulence onto the ion scale is minimal 
Ref. ( IWCF07( IGJOSbp .) The turbulent frequency is defined u = Vthi/L relative 
to the ion thermal speed. It is important to make these assumptions explicit at 
this point because they will be implicit during the analysis. In the case of electron 
dynamics (for ion-scale turbulence), we will be able to exploit an additional small 
parameter: the electron-ion mass ratio is taken to be the order of the expansion 
parameter: 



Formally, this will require fractional ordering of the equations when square roots 
of the mass ratio appear. However, there is nothing fundamentally new about 
the analysis and the results follow easily from the ion case, as we will see. 

2.4 Averaging and scale separation 

As previously mentioned, the multiple scale approach aims to obtain both tur- 
bulent and classical (or neoclassical) theory simultaneously. To this end, we use 
smoothing averaging to separate the equations for equilibrium-scale quantities 
from the rapid turbulent dynamics described by gyro-kinetic theory. This section 
introduces time and space smoothing operators which averaging over regions of 
space and time intermediate to the cyclotron, turbulent and transport scales. It 
should be stressed that these are tools of the formalism and the final equations 




(2.9) 
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do not depend on a specific 

2.4.1 Gyro-average and single peirticle motion 

The gyro-average exploits the disparity between the cyclotron frequency and the 
turbulent frequency The rapid rotation of individual particles sample turbulent 
electromagnetic fluctuations which are static on this timescale. This means that 
turbulent timescale dynamics are influenced by the orbit or gyro-average of the 
fluctuating flelds. Let's examines individual particle motion more closely The 
velocity vector is defined using cylindrical coordinates aligned with the magnetic 
field. 

V = I'ljbo -I- i'_L(cos'j?ei -I- sin'j?e2). (2-10) 

where the unit vectors bo, ei and e2 form a local right handed coordinate basis 
i.e. ei X e2 = bo- (Because these basis vectors are defined by the equilibrium field 
geometry, it should be apparent that they vary on the macroscopic, L, spacial 
scale and the slow, r, time scale.) The fastest motion is the gyro-motion: 

M 

— ^-Qo + O{eno). (2.11) 

where Qq — QBo/m is the cyclotron frequency for the species of interest. An 
individual particle's velocity is thus fiuctuating rapidly on the fast time scale. If 
we look at the motion of the gyro-center, however, one finds it's fluctuations are 
small. One can use the gyro-average to separate smooth turbulent scale motion of 
the gyro-center to obtain the drift behavior of particle motion in the gyro-kinetic 
limit. The gyro-center position is defined by the the Catto Transformation 
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R = r — p 

where p is the larmor radius vector defined by 



(2.12) 



bo X V 

fin 



(2.13) 



The gyro-average is defined as follows 



(A(r,v,t))^ = i-^ '^A(R- ^^,v,t)rf^. (2.14) 

It is an average over the gyro-angle d with the gyro-center position R held fixed. 
Note that the parallel velocity, magnitude of the perpendicular velocity and all 
equilibrium quantities are fixed since they do not significantly vary over a parti- 
cle's trajectory during a gyro-period. If one applies the gyro-average to the rate 
of change of the gyro-center position, dH/dt, one can obtain the drift velocity of 



a particle's gyro-center (see section 2.6.3). To reiterate, this quantity varies on 



the medium turbulent time scale, with cyclotron time scale behavior having been 
smoothed away. 

2.4.2 Spatial averages 

To separate turbulent spatial scale from the equilibrium scale we employ a patch 
volume average. The patch volume Vp is defined around a point (x, y, z) defined 



using locally fiat coordinates as indicated in Fig. 2.1 



AxAyAz ^^-^^^ 
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Figure 2.1: Patch Volume and Scale Separation: The patch volume Vp (colored 
brown) is chosen to be small enough such that the equilibrium geometry is fiat 
but large enough to contain many turbulent scale (correlation) lengths. 

where the intervals Ax, Ay and Az define a region with dimensions intermediate 
to the microscopic and macroscopic spatial scales-i.e. we require 



p<Ax, Ay, Az<^L. (2.16) 

Formally, we take Aa; ~ Ay ~ Az ~ \/eL. The details of the patch average 
definition are less important than its practical utility. As stated, it smooths 
away turbulent scale dependence. For example, if one assumes that |VA| ~ 
0{A/ p) it is easily shown that |V(P[A])| ^ |Vv4|. In analogy to the gyro- 
average which smoothes out fast time scale behavior, the patch average smoothes 
microscopic turbulent spatial structure. We will use it often to separate quantities 



into equilibrium and turbulent parts (see Sec. 2.5). 
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2.4.3 Time average 



The last average acts between the slow transport time scale and the medium 
turbulent time scale, separating out turbulent time scale structure. The time 
average is defined 

l't+At/2 

A = — A dt', uj-'^ < At < T (2.17) 

Jt-At/2 

Formally, we take At ~ e^^^r. 



2.5 Maxwell's equations and potentials 

In this section we demonstrate use of the ordering assumptions to obtain some 
simple results from Maxwell's Equations. Let's start with Faraday's Law 

= -V X 5E. (2.18) 
dt ^ ^ 

Upon ordering, the dominant equation is 



V X 5E = (2.19) 

This means that the inductive part of the electric field must be higher order than 
the electrostatic part so we may write 

-W{ip) + 0{e^) (2.20) 
This suggests using potentials to express the fields. Thus we define 
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A = Ao + (5A 



(2.21) 



where 



Bo = V X Ao (2.22) 
5B = V X 5A (2.23) 

For definiteness we use the coulomb gauge, V • A = 0. Note that our ordering 
requires 5 A ~ 0{e^Ao). Additionally, it should be apparent that Aq is an 
equilibrium quantity and thus has only slow scale dependence like Bq. Likewise, 
6 A has turbulent scale dependence. The perturbative electric field becomes 

dA 

SE = -V(^ - — (2.24) 

where dA/dt is composed of equilibrium and perturbative parts of the same order 
in size. We now show that SE and 5B both have small patch averages. Since SE 
is mostly electrostatic we may write 

V[SE] - -V[Vip] (2.25) 
To demonstrate use of the patch average, we apply it to the electrostatic field: 



recaUing that p Ax, Ay, Az <^ L. Then we have 



p[iE]~0(-x|--y^-z£)«^E (2.27) 
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The proof works likewise for 5B. 



2.6 Ordered Fokker-Planck equation for ion species 



In Sees. |2.6] and |2.6.4] we will be solving the ordered Fokker-Planck equation with 
the goal of arriving at the gyro-kinetic equation and classical equation which are 
used to solve for 6fi. The ion case is treated in more detail because the electron 

we 



case can be approached analogously. Using the orderings listed in Sec. |2.3 
write down the Fokker-Planck equation with terms ordered in powers of e relative 
to VthiFo/L. Each set of terms with the same power of e will yield an equation 
to be solved. In this section, we will solve the Fokker-Planck equation for ions at 
each power of e for the three largest orders. 



^ + v-V/+^(E + vxB).§^ = /) (2.28) 
ot m av 



2.6.1 0{e-^) equation: 

One term in the Fokker-Planck equation is larger than all others. 

from which we deduce that Fq is independent of gyro- angle i? so that: 

Fo = Foir,v,v^,t) (2.30) 

Now we proceed to C(l): 
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2.6.2 C(l) equation: 



The next order ion equation is 

V • VFo + ■ V5/i + ^ (-Vy. + v x 5B) ■ ^ - (^] = C{Fo, Fq) 

(2.31) 

Applying the patch average to this equation, we obtain a separate equation con- 
taining the slowly terms. 

V ■ VFo + -Qo (^^|y^) = C{Fo, Fo) (2.32) 
From this equation we use Boltzman's H-theorem to show that Fq must be 



Maxwellian. To do this, we multiply Eqn (2.31) by 1 + logFo, integrate over 
velocity space. We find that upon integrating, most of the terms on the left side 
vanish, leaving us with 

J ciVi;||bo ■ VFo = J d\C{Fo, Fq) IuFq = (2.33) 

We make the additional assumption that bg ■ VFq = 0. This assumption is valid 
for systems with closed fiux-surface geometry (as proven in the Appendix) such as 
toroidal fusion devices because rapid transport along the equilibrium magnetic 
field lines enforces constant equilibrium temperature and density on each fiux 



surface. What we are left with from Eqn (2.33), to 0{1), is 



d\{\nFo)C{Fo,Fo) = (2.34) 



Eqn (2.34) and the H theorem tell us there is no entropy production at this 
order. Our first significant result in this section is found because we know that 
the entropy production can only be zero if Fq is a local Maxwellian. 
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n 



exp 



(2.35) 



Given this result, we can obtain information about 6fi by using our solution for 



Fo in Eqn (2.31). 



■ V5/i - no (^) = -v^ ■ V(^) Fo - • VFo (2.36) 

We have dropped the factor v^\hQ ■ V(^) Fq from the right hand side because it 
is negligible at this order. The following particular solution satisfies the eqaution 
at this order: 



5/ip = -(^)Fo-p-VFo. (2.37) 



where p is the larmor radius vector defined by Eqn. 2.13 The first term is the 
perturbed Boltzmann response of the particles to the potential. The potential is 
static on the gyration time scale and therefore the energy S = {l/2)mv'^ + q(f is 
conserved to this order. The second term can be recognized as the beginning of 
an expansion of Fq about the gyro-center position: 

Fo(R) ~Fo(r)-p- VFo(r) (2.38) 
We absorb these terms into a new definition of Fq in gyro-center coordinates, i.e. 



Fo = F^(R)exp-(^) (2.39) 
where F^ is the maxwellian distribution with real-space coordinates 
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n{t,r] 



3/2 



exp 



(l/2)mf2 

nt,r) 



(2.40) 



To obtain the homogeneous solution for Sfi in Eqn. (2.36), we rewrite the left 
hand side of this equation using the following identity: 



^^0 



R 



■ V 



(2.41) 



Dropping the right hand side of Eqn. (2.36) will give us the equation for the 
homogeneous part of 6fi, Sfih- 



dSf 



Ih; 



0. (2.42) 



R 

Thus the homogeneous part is independent of gyro-angle at fixed R i.e., 



Sfih = h{R,v,v^,t) (2.43) 

Sometimes h(R,v,v±,t) is called the Guiding center distribution. As we 
show in next order h satisfies the gyro-kinetic equation. Combining all 0{1) 
results gives the following distribution function: 



/(r, V, t) = Fo(r, v, t) + h{R, v, v^, t) + 5h{v, v, t) + ... (2.44) 

This fixes the form of the solution for the distribution function /. However, we 
still need to derive equations to describe the evolution of /;.(R, f , f_L, t), n(t,ip) 
and T{t,ip). Now we proceed to 0{e) where we obtain the gyro-kinetic equation 
and classical equation to solve for h: 
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2.6.3 0{e): Gyrokinetic equation 



It is now straiglitforward to take Eqn (2.44) and Eqn (2.40) and proceed to 



the next order in the Fokker-Planck equation. We would then gyro-average this 
equation to obtain the gyrokinetic equation. This process is somewhat tedious. 
We find an improved route by recasting the Fokker-Planck equation in gyro-center 
variables: 



V X bo ^ 1 r, mv 
R = rH — , c = i^mv + qip, fi 



2Ba 



The Fokker-Planck equation takes the form 



dl dR dl_ dfxdl d£dl d^dl_ 

dt^ dt ' dR^ dt 9/i ^ dtdS^ dt dd ~ ^ ■ ' 



Our final equation will involve the gyro-average defined in Eqn (2.14) of dR/dt 



dS/dt, dd/dt and d^/dt. Calculated to necessary order, these quantities are 



<f >R = ^l|bo + W X + ^bo ■ Vbo + f VlnBo) (2.46) 

<f>K = ^(^-vr^) (2.47) 

(f>R = -^o (2.48) 

(f>K = (2.49) 

where we define x = — v ■ 6A. To be explicit, we note that in these variables 
/ = -Fo(R, £, t) + h{R, £, /i, t) + 5/2(R, £, fi, ^, t) + ... The (9(e) equation is found 
to be 
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dh I dR d ( v I T? \ i dn dh 



dt 



dt d^i 



C{h- p-WF^.F^) 



(2.50) 



inside the collision operator refers to Eqn. 2.35, and we will write C{.,Fm) 
as C{.) from now on. Now we can eliminate 5/2 to obtain an equation in h only 
- to do this we gyro-average the equation. With some algebra we arrive at the 
following. 



VTo I dt 



(vd + V J ■ VFo 



(2.51) 



where the equilibrium and perturbative parts of the drift velocity are split as 
follows 



X UJbo-Vbo + ^VlnSo 



Bo 



X (Vx)r 



(2.52) 



Note that the term C{—p ■ VFm, Fm) in Eqn. 2.51 was annihilated by the gyro- 
average. Hidden in this equation are actually two independent equations. This 
is due to the multi-scale dependence of h. We recall the patch operator and time 



average from Sec. 2.4.2 and define their action on h as follows. 



V[h] = ha 



(2.53) 



This allows us to split up h into two parts: h = hd + hgk, where "cl" stands for 
classical (which will be the neoclassical response in toroidal geometry) and "gk" 
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stand for gyro-kinetic. Applying V and the time average to equation (2.51) gives 
the equation for the "classical" part of h (an equation of slow varying quantities). 



v\\hQ ■ Vhci - C{hc 



dt 



(2.54) 



Now subtracting equation (2.54) from equation (2.51) above gives the Gyro- 
kinetic equation: 



^ + (^;|ibo + + ■ Vh,, - {C{h,,))^ = - V, ■ VFo (2.55) 

In some loose sense the Gyro-kinetic equation is the kinetic equation for rings 
of charge centered at R(it) of radius v±/fl. 



2.6.4 Kinetic and classical equations for electron species 



As described in Sec. 2.3.1, we assume gradients and time derivatives act on the 
electron distribution function at the same order as ions. The significant difference 
when deriving the kinetic equation for electrons is the smaller mass. When re- 



examining terms in the Fokker-Planck equation, Eqn (2.28), we order the square 



root of the mass ratio as ymg/mj ~ ^/e which formally introduces fractional 
ordered terms in our perturbative expansion. However it is much more convenient 
to take a maximal ordering approach and include root mass ratio terms until the 
final equation (the kinetic equation for electrons) is obtained at which point we 
can introduce the effects of the mass ratio as a subsidiary ordering. Thus we 
can immediately write the gyrokinetic and classical equations from the previous 
section for the electron species and include only dominant terms in the root mass 
ratio: 
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v\iho ■ Vhci - C{hci) = -?^V|| • (2.56) 

Iq at 

and 



F dSA. 

(v\iho + Ve) • Vhgk - C(hgk) = -Ve • VFq - g-^V|| • — (2.57) 

io ot 

where Ve is the perturbative magnetic field drift velocity for electrons 



60 X V(v\\SA\\) 
Ve = ^ — ^ (2.58) 

-Do 



2.7 Maxwell's equations for gyrokinetics 

The gyro-kinetic equation from the previous section determines how to evolve 
the turbulent distribution function. This section describes the equations that 
will evolve the turbulent potentials, and consequently, the turbulent fields. 

2.7.1 Queisi- neutrality 

+ 27rg/ / v^dv^dvw {hi,gk{R,v^,v\\,t))^ 
= ^ + 27re// v^dv^dv\i {he,gk(R,v^,v\i,t))^ (2.59) 



2.7.2 Parallel Ampere's law 

bo • A = /xo J|| = J d^vvii [q fi-e fe 
which results in two equations upon patch averaging. 



(2.60) 




V 6A\i = /xo27r / / dv±dv\iv\\v±[q hi^gk - e he,gk] (2.61) 
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and 




bo • V^Ao = /io27r / j dv^dv\\v\\v^[q hi^a - e K,ci\ (2.62) 



2.7.3 Perpendicular Ampere's law 

(V X B) X bo = /^oy c?'v(v X bo)[g fi - e Q (2.63) 
again, giving two equations. 

- Vi_5B\\ ^1^0 J c?^v(v X ho)[q hi^gk - e K,gk\ (2.64) 

and 



5obo-Vbo--V^52=/io j dMvxBo)[g(-p,-VF„0-e(-Pe-VF„e)] (2.65) 
2.7.4 Pressure balance 

The right hand side of this equation, /ioJo x Bq, can be manipulated to give 
the familiar equilibrium force balance. The surviving terms after velocity space 
integration can be identified as the polarization current (recall that we have 
assumed there is no equilibrium scale electrostatic potential). 
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Jo X Bo = 5o J rf3v(v X bo)[g(-p, ■ VF^) - e(-p, ■ VF^e)] 

= j rf^v(v X bo)(v X bo) ■ V[miEra^ + nieFrae] 

f 

= / d\^V±[miFmi+meFme\ 

= VPo (2.66) 



where we recall = and used the fixed-r ring average between lines one and 
two and use that -Pq = | Jl'^^Fmi + ^^-F^e]- Substituting this into Eqn. 



2.65 



we arrive at the following equation 

Eo'bo ■ Vbo - l^±B^o = /xoVPo (2.67) 

This equation describes the macroscopic equilibrium of the plasma. This com- 
pletes the solution on the medium timescale fluctuations and we now go one more 
order to obtain the slow evolution of the equilibrium. 



2.8 0{e^): Transport equations in the slab limit 

To proceed to the problem of transport and evolution of the equilibrium, a full de- 
scription of the equilibrium geometry is necessary. The full axisymmetric toroidal 
case is treated in the Appendix. Here we calculate transport in an illustrative and 
simple limit of gyrokinetics, the straight field line "slab" geometry. The equilib- 
rium magnetic field lines have a constant direction {z), and the curvature drift is 
eliminated from the gyrokinetic equation. The x coordinate is taken as the flux 
surface label so the equilibrium has variation only in this direction. Also note 
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that for this case the classical part hd will be neglected as will the the evolution 
of the background magnetic field (dAo/dt). The system is taken to be periodic 
in the y and z directions. In place of the time average and patch average we will 
use a single smoothing operator which will combine an intermediate-scale average 
in the x-direction with a full average over the system in the y and z directions 



and the intermediate time-scale average defined earlier by Eqn. 2.17 



where the system extends a distance L in the y and z directions and, thus, 
V = AxL"^ is the volume over which the smoothing operator averages. To obtain 
the slow time dependence of uq and Tq we consider the moment equations of 
the full kinetic equation and apply the smoothing average. Many terms will be 
zero outright and we can avoid tabulating all the 0{e^) terms from the expanded 
equations. Also, since the moments are taken in real-space coordinates, it will 
often be convenient to refer to terms explicitly in these real-space coordinates 
(as opposed to gyro-center coordinates). I.e. we will make reference to 5fi 
which we recall is the first order correction to Fq in non gyro-center coordinates: 
/ = F^(r, v) + 6h + 0{t^) where 6h = -p ■ VF^ - q^Ejn + h. 

2.8.1 Particle transport 

To obtain the time evolution of the density, we will integrate the Fokker-Planck 
equation over velocity and apply the smoothing operator . 

S[ J d\^ + V . V/. + ^(E + V X B) ■ = d\Cif, /)] (2.69) 
The integral of the collisions over velocity will be zero to conserve particles. The 
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(E + V X B) • terms will be zero because they are a perfect divergence in 
velocity space. What is left is the continuity equation. 

S[Jd\^ + V-{^rf,)] = (2.70) 

The y-derivative and z-derivative parts of V ■ (v/^) will spatially average to zero 
due to periodic boundary conditions in those directions. Finally, the x component 
of the V • (v/s) term can be rewritten by considering the following: 



/ d\{v X bo) • X— ■ (v X bo/s) 

= - y d^v [(v X bo) X bo] • xfs 
= J dM^-x)fs 

Thus, we can write 



f d 

/ d^v{v X hofs) ■ ^(v X bo) • X 

(2.71) 



d_ 
dx 



dM'Vxfs) 



f d 
5x • / rf'V(v X bo)^ • (v X bo/, 

J rf^v(v X bo 



(v X bo) 



dv 



— (v X Bo) • ^ 
m ov 



(2.72) 



where we used ^ • (v x bo) = and introducted the notation d:x. — xd/dx. We 
can rewrite [^(v x Bq) ■ ^] using the full Fokker-Planck equation. 
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5. 



fv X bo 



3,.(v X bo) 



(v X Bo) 
dt 



vV/. + ^(E + vx5B)-|^ 



Cif, f) 

(2.73) 

The advantage of writing the v- V/^ term in this form is that the factor of '■'^^^"^'^ 
in front of the Fokker-Planck terms will reduce the order of the terms inside since 
{Q^^ ~ 0{-)). This is as far as we can get before using the solution for fg. Now, 



we can rewrite Eqn 2.70 as 



1 

5x ■ / rf^v(^5[^ + v_V^+ ^(E + V X 5B) . ^-C{f, /)] (2.74) 

"-^v^ 3 ^ ^ ' 

2 4 



If we plug our solution of fs into term 1 of Eqn 2.74, we obtain up to 0{e 



We identify J d^v^^^ = as the quantity we are solving for, namely the time 
derivative of the equilibrium density. The 5fi terms are negligible as follows: 
firstly, the boltzmann response vanishes after space averaging (we have assumed 
there is no equilibrium-scale electrostatic potential) and second, gyro-center dis- 
tribution function h contributes negligibly after smoothing (note that we have 
neglected to include in this derivation but contribution from this term would 
be obvioiusly be an order epsilon smaller than dFo/dt). Finally, the 5/2 term in 
Eqn 2.75 will drop out after performing the time average. Term 2 in Eqn |2.74 is 
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treated analogously to term 1 but the order of each term is one smaller because of 
the preceding factor of ^'^^q"'*'^ - Thus, it contributes negligibly. The third term in 



Eqn 2.74 has derivatives in y and z which vanish by periodicity after averaging. 



What remains is the x-gradient, 

S[ I d^v ^^^y •^ t;.^)(Fo + h + Sm (2.76) 

None of these terms will enter at 0(e^). To show this we first note that h and 
Fo, are functions of gyro-center position R. These terms can be first spatially 
averaged over y. This removes any odd dependence in due to gyro-center 
dependence. The resulting terms are odd in Vx and can be integrated away. The 
5/2 term will drop in ordering after averaging it over space. The fourth term in 



Eqn 2.74 can be written to 0{e ) as 



I ^3^ (vxbo)-£ ^(E + V X 5B) . ^(Fo + h)] (2.77) 
J il m av 

The parts of Fq are negligible as follows: The Boltzmann response term is zero 
under velocity differentiation. The gyor-center correction —p ■ VF^ is an order 
smaller after smoothing over the turbulent fields E + v x (5B. The Maxwellian 
part Fm is argued away by a combination of periodicity in y and z, time averaging 
and oddness in velocity velocity space. We now write E + v x (5B as — Vx — 
To simplify the h term, we integrate by parts in cartesian velocity coordinates, 
so that 



S[ [ d^'^^i-Vx - ^) ■ ^h] =S[! d\^x ■ (Vx X bo)/.] (2.78) 
J \l m at o\ J Bo 

Where the ^ term is removed by by noting that ^ is small in our ordering and 
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employing Eqn 2.41 after integration by parts to remove the rest. The final term 
in Eqn 2.74 simply involves plugging in values for the distribution function. 



(2.79) 



Recall that the Maxwellian is notated so as to distinguish it from the modified 
Maxwellian Fq defined in Eqn 



2.40 



Collecting terms, and recalling the in front 
of the integral from Eqns 2.78 and 2.79 from Eqn |2.74 , we have the transport 
equation for particles: 



dt 



-Do 



(v X bp 



-C{p-VF^,F^) 



(2.80) 



2.8.2 Heat transport equation 

To calculate heat transport , we multiply the Fokker-Planck equation by ^mv^, 
integrate over velocity, and smooth. Following the same methodology as used 
with the particle transport. 



2 at m o\ 



mv 



d^^ S[CifJ)] (2.8i; 



The first term of Eqn 2.81 yields 



3d(T„n„f 



2 m +5|/<;vl™v=|(*AH-v,)i 



(2.82) 



The first term is equilibrium pressure evolution fTono = / d^v'^F^. The other 
terms in the time derivative will be negligible after smoothing. The second term 



in Eqn 2.81 can be rewritten in the same manner as Eqn 2.71 
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dx / 2 



5x 



rf=^Vy(v X bo)^ ■ (v X hofsV^) 



fv X b 



0) 



n 



— (v X Bo) ■ {v — 
m av 



(2.83) 



The last line was obtained using (v x bo) ■ ^ = 0. Once again, we are able to 
substitute in the Fokker-Planck equation, and obtain 



m 



fv X b, 



oj 



n 



^(v X B„) . 

m av 





2 



3 2 (V X bo) 

f 1^ 



^+vV/. + ^(E + vx5B).|^ 



C{f,f) 
(2.84) 



The difference between this equation and Eqn 2.73 is obviously the f ^ inside the 
integral. Fortunately, all the terms that were negligible in density transport will 
be negligible in this well, because the arguments relied on parity in 

and Vy, spatial smoothing and ordering, and integrations. Thus, what we really 
need to calculate is 



51 / #v !^|^(-Vx 

J 2 \l m 



dA dh. 



(2.85) 



2 Vt m dt dw 

Here, the integration by parts will be a little messier than with the analogous 
term in the particle transport equation. The ^ term is negligible as before. The 
result is 



9. 



,o ^r"^f^(Vy X b) , m(vxbo),, 

'-h + ^((v ■ V)0)/i] 



2B. 



(2.86) 
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We will later see that the second term of Eqn |2.86 will cancel. The collisional 
terms of Eqn |2.84 work similarly as with density transport. 



-d^-J rfW^^^^^^^C(-p ■ VF^, F^) (2.87) 
Unlike the case of particle transport, the electric field term will not integrate 



away from Eqn 2.81 First we manipulate as follows: 



3 mv^ q df 



The scalar potential part of the electric field is an order larger than the ^ term, 



but as we will see most of it will cancel. This will result in the part of the scalar 

3A 

at ■ 



potential that does not average away acting on the same order as 



/ rf^v [gv ■ V0/] = - / rf^v [gf /] + / d^^r [g(f + v ■ V0)/] 
= - / ci^v [gf /] - / ci^v [q{% + v ■ V/)0] + / rfV [g^^ 

+ / d^YqV{(j)f) ■ V 
/ rf^v [gg/] + g / rf^v . [(E + V X B)0/] - / rfV [g(C(/, /))0] 



y at 



+ /d3vg(vV)(0/) 

+ /d3vg(vV) (</)/) 



(2.89) 



We have integrated by parts in time and space between lines one and two. Be- 
tween lines two and three we substitute the full Fokker-Planck equation ^ = 
C{f,f) and, finally, we use gauss's law in velocity space and that collisions con- 
serve particles to eliminate the velocity divergence and collisional integral to 
obtain line four. The second term on the final line will be negligible after time 
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average. We find that the very last term will cancel the last term of Eqn 2.86 
when we compute the contribution of h from / in the equation. 



I' 



d''vq{vV){(t>h) 



d_ 
dx 



(f\qv^{(j)h) 



d f ,3 ./V-V, fdh\ . 



dx 

J -Do 

Between lines 1 and 2 we integrated by parts in velocity coordinate 19. Between 
lines 3 and 4 we have integrated by parts in space to move the gradient onto the 
(f). We now want to combine our results from the electrostatic part of the electric 



field in Eqn 2.90 with the inductive part. 



m 
~2 



where the Maxwellian part of Fq averages away due to velocity parity. We now 
plug in our solution for 6fi. The adiabatic/Boltzman term gives 



S[ I d'^rq^^[<P - - . AJ(-g|F^)] = (2.92) 

since: the A\\ term is odd in f y and the A_i_ gyro-averages to zero (or equivalently 
is odd in v±). What is left will time average away, leading to the result that at 
this order the E ■ v^F^ term produces no heating. Returning to the original 
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moment equation in this section, Eqn. 2.81 the v x B terms will not contribute, 
as magnetic fields do no work (it can be written as a perfect divergence and 
eliminated by integration over velocity space). The final term we need to calculate 
is the coUisional energy exchange terms between species are now included - note 
like particle collisions do not produce a loss of energy and thus do not appear 
(this is easy to prove). The collisional energy exchange is standard since to this 
order it is between Maxwellian species. Specifically: 



J d\]^mv^C,,{fi, /e) = n,vf{T, - T,). (2.93) 

where is found in ( ]HS02[) - it is ^/mjrrii smaller than the ion-ion collision 
rate which is itself \fm^rni smaller than the electron-ion collision rate. We now 
collect all the results and write down what we have so far for the heat transport 
equation. 



3 9nojT( 



Oi 



2 dt 

/• 3 mt;2/ (Vxxbo) (v x bo) A 
■ / V— \ S\ h\ + ^ [C(p • VF™, F^)] 1 

+ y" S^qS^-^h\ + n^v\\T^, - To,) (2.94) 

We can rewrite the gyrokinetic heating term qhdx/dt (this is a necessary step for 
establishing entropy balance later). We write the smoothing operator explicitly as 
a spatial average (ignoring the time average for now) and use Parseval's identity 
to write the integral as a sum of Fourier components so that the gyro-angle 
average may be shifted to x- The result is 
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r^jdw-^h 



(2.95) 



Now using the gyrokinetic equation, Eqn. 2.55 this expression becomes 



(Vx)RXbo 
Bo 



VF,f-j^^jd^.f{C{KF^))^ 



dt V 2Fo . 



(2.96) 



where we have multiphed the gyrokinetic equation by and apphed the indi- 
cated averages. The final term will average away upon time-averaging. We may 
reverse the change of variables trick and revert back to the smoothing operator 
form to plug into the heating equation. 



Thus, combining all the surviving terms from Eqns 2.82, 2.86, 2.87, 2.96, and 



2.93 will give us the heating equation. 



3 driQiTQi 
2 dt 



■ J d\ i^S[ '-^h] + \^ " [C{p ■ VF„, F^)] 



,m(\ X bn 



2Bn 



2n 



[ d\S[^^^ ■ VF^^h] - I d\S[^C{h,Fj] 



+ no4'{T,-T,). 



(2.97) 



58 



2.8.3 Entropy balance 

Establishing entropy balance is a good way to ground the results of the gyrokinetic 
transport calculation to an exact result from statistical physics. A tokamak (or 
any other plasma) is ultimately a coUisional (dissipative) system. The generation 
of entropy by the collision operator provides the mechanism for the ultimate 
dissipation of turbulent fluctuations and in particular prevents the formation of 
unphysically small scale features of the distribution function. Also, calculation of 
entropy balance has the added benefit of a consistency check to confirm that the 
terms in our transport equations have been calculated correctly. The time rate 
of change of the mean entropy per unit volume (contained in the volume over 
which we are smoothing) is given by 

^ = -yrfV5[^/l./l (2.98) 

We multiply the Fokker-Planck equation by 1 + In / and integrate to obtain the 
relation 

J d\^-^ + j d\S[^-V{f In /)] = j d\S[C{f, f) In /] (2.99) 

We will now use heat and particle transport equations to demonstrate that en- 
tropy balance does indeed take the form 

^ = 9. • y" d\S[Yf\^f] - j d\S[C{f, /) In/] (2.100) 
It is not hard to show that ^ is mostly due to the background Maxwellian, i.e. 

~dt^~dtj ^2.101) 
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Performing the integrations give 



as 

~dt 



3 TTli 

Inrin H — ln( — — ' 



drio ^ 3 9noTo 



dt 



2T dt 



(2.102) 



We can now substitute in from Eqns 2.80 and 2.97 



as 
at 



- J rf^v (in no + I In(^) + l) 9. ■ [s[^^h] -pC{p- F^, F„ 
0^ . J d^v^ (S[^h] +pC{p- VF^, F„ 



■},[jd^^S[^h].V{\nF^)T + jd^^S[f;^C{KF^)] 



(2.103) 



What is left to do is some integration by parts and arrangement of terms followed 
by evaluation of VlnF^, on the third line. The result is 



f = 9x • / rf3v(ln F^ + l)(5[Mo /i] -pC{p- VF^, F^)) 
- / dMP ■ V In FjCip ■ VF^, FJ - J d^wS[^C{h, F^)] 

+noz/i"^^ (2.104) 

which completes entropy balance. Notice that the finally two lines can be easily 
seen to come from j d^\ In f C {f , f) . The two terms on the second line repre- 
sent the nonvanishing parts of J d^\^pr^C{6fi, Fm) and the final term is simply 
/ d'v\nF^CiFm,i ) Fm,e) 

2.9 Poynting's theorem 

A standard derivation of Poynting's Theorem begins taking the inner product of 
Ampere's law with the electric field 
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d'E 1 

J-E = -eo^ ■£+ — E- (V X B) (2.105) 

One continues by applying the identity V-(AxB) = B- VxA — A-VxB 
followed by substituting Faraday's Law V x E = — ^ to obtain the standard 
form of Poynting's Theorem: 

J.E=-l(l^+.„f!)-lv.(ExB) (2.106) 

2 V/io ot dt J 

This gives a complete picture of the electromagnetic energy flow, making no 
ordering assumptions, but throws away details by lumping together all flelds 
(electrostatic, inductive, equilibrium, perturbative). We can obtain Poynting's 
Theorem in three separate parts and use our ordering assumptions to get a de- 
tailed picture of electromagnetic energy flows. We do this by taking the inner 
product of ampere's law separately with three different parts of the electric fleld: 
the electrostatic part (Eg = —V(p), the equilibrium inductive part (Eio = — ^^) 
and the perturbative inductive part (Eii = — ^^)- We proceed exactly as with 
the standard Poynting's Theorem derivation above. Only lowest order terms 
are retained in each equation. We order the J ■ E terms using preceding kinetic 



analysis (see Eqn. A. 73) and flnd they are each 0{e uT). We then use a patch 
volume average and time average to separate the multi-scale dependence. First 
the electrostatic part we flnd 

J-E, = V- (E, X B) (2.107) 

Where we have neglected the displacement current term which will enter at a 
higher order in Poynting's Theorem for both non-relativistic and e ordering ar- 
guments. We have also used here that V x Eg = because Eg is the gradient of 
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a scalar. Wc note that J ■ is an order smaller than the remaining term on the 
right hand side of the equation. This gives us the following result, at dominant 
order: 



• (V X B) = (2.108) 

We can interpret this to mean that the electrostatic Poynting flux is everywhere 
zero. At next order we learn something about J • E^: 

J-E^ = -— V-(E^x5B + E,2xBo) (2.109) 

where E52 is the second order electrostatic potential - for which we do not solve. 
However, all terms on the right hand side are linear in the fluctuating fields, 
so applying the patch average, they are eliminated (increased in order) and we 
obtain the result 



V[J • E^] = (2.110) 

which means that the electrostatic field does no work on the total current. Now, 
the equilibrium inductive part: 

1 1 dB'^ 

J • E,o = V • (E,o X B) - --° (2.111) 

Ato 2/^0 ot 

Again, neglecting the displacement current and higher order terms. We apply 

the patch average to remove the 6B perturbative term on the right hand side. 



1 1 dB^ 
J • E,o = V • (E,o X Bo) - ^ (2.112) 
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This equation describes the balance of energy flow due to Ohmic heating. Lastly, 
the perturbative inductive part of the electric field. 



3-Ea = --V ■ (E,,i X B) - -Bo ■ ^ (2.113) 

The patch average and time average of the right hand side removes all terms due 
to ordering arguments. 



V[3 ■ Ea] = 



(2.114) 



Thus from this result and Eqn. |2.110 we find that the fluctuating electric field 
(both inductive and electrostatic) does no mean work on the total current. 
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CHAPTER 3 



Primary and Secondary Mode Theory 

The following statement applies to the figures included in this chapter, which 
have been taken from the paper \Plu01 ): Reprinted with permission from Gabriel 



Plunk, Physics of Plasmas, Vol. I4, Issue 11, Page 112308, 2007. Copyright 
2001, American Institute of Physics. 



3.1 Introduction 

A natural first step in understanding anomalous transport is the calculation of the 
primary instabilities which drive micro-turbulence. We will begin this chapter by 
reviewing a fluid description of the relevant instabilities driving core turbulence 
and also solve the gyrokinetic linear dispersion relation. This chapter continues 
with secondary instability, a powerful and physically intuitive analytic theory that 
helps to explain mode saturation and transition to turbulence. The material on 
secondary instabilities is copied to a large extent from the author's work (lPlu07[) 
with permission from the publishers. 

The ion temperature gradient (ITG) drives modes at spatial scales compara- 
ble to the ion Larmor radius and is typically the strongest source of transport. 
Electron thermal transport due to turbulence driven by the electron temperature 
gradient (ETG) and/or trapped electron modes (TEM) is regarded as a strong 
second candidate. This chapter investigates ITG and ETG turbulence. 
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In this chapter and the following chapter we take the simplification that the 
electrostatic perturbations are dominant over the magnetic field perturbuations. 
In this chapter we also assume that the coUisional terms are subdominant. In 
plasma turbulence, collisions ultimately must play a critical role even when they 
are very weak and do not enter at the scales of the primary instability. As 
with fiuid turbulence, the steady state of a turbulent plasma is characterized 
by injection and ultimate dissipation, without which entropy production and 
velocity space distribution would develop unphysical characteristics. This issue 
will be explored in more detail in the following chapter. 

The system is closed by the quasi-neutrality condition and, for most of the 
chapter, the Boltzmann or "adiabatic" response will be assumed as a simplifica- 
tion. In particular this means that for ETG driven turbulence, ions are assumed 
to have a Boltzmann response and for ITG driven turbulence, the electrons have 
the appropriate Boltzmann response. The previous chapter has built the founda- 
tion for a local treatment, with spatial derivatives of equilibrium quantities taken 
as constants. As demonstrated later, the local approximation is remarkably accu- 
rate in calculating linear growth rates of micro-instabilities. In nonlinear regime, 
we expect the local case to give insight into tokamak turbulence and provide the 
basis for a more complete treatment in future work. 



3.2 Normalized equations 

It is convenient to now use a normalized version of the gyrokinetic equation. 



Eqn. 2.55, A locally fiat orthogonal coordinate system is used where, as is the 
convention, x is identified with the radial direction, z is the coordinate along the 
magnetic field and y is the binormal coordinate (normal to both z and x). The 
gradient lengths are taken to be constants and the magnetic field is taken to have 
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constant curvature (i.e. b- Vb and Vlni? are constant vectors where they appear 
in the drift velocity below). The normalized gyrokinetic equation is written 



dh dh 



(3.1) 



with 



K = x{ri~^ + v^/2- 3/2) 
vd = zx (vlLrh ■ Vb + V In B) (3.2) 



with T] = Ln/Lx- The bracket average in Eqn (3.1) denotes the gyro-average, 
an average over gyro angle at fixed gyro center R. In this normalization the 
gyrocenter variable is defined H = r — p = zz + Yy + Xx = zz + {y — v± sin '&)y + 
{x + v± cos'&)x. The normalized quasi-neutrality condition is 




27r / / v±dv±dvii (h)^ = (1 + r)Lp - {5{itg)) t'^ (3.3) 

with 



The fiux surface average ip = f f ip dydz / {LyLz) is an average over the periodic 
box with sides Ly and Lz. The notation {.)^ is used for an average over gyro- 
angle ^? at fixed spatial coordinate r. Note that the fiux surface averaged potential 
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enters crucially for the Boltzmann (adiabatic) electron response in ITG turbu- 
lence but is absent for the Boltzmann ion response for ETG turbulence. This is 
the only difference between the equations governing electrostatic ITG and ETG 
turbulence. This results from the different physical mechanisms that create the 
boltzmann response. For ETG, the ion Larmor motion provides rapid motion 
perpendicular to the field lines that allow the ions to "sense" electrostatic per- 
turbations. In the case of ITG, the Larmor radius of electrons is small compared 
to the fluctuations and it is the rapid motion along the fleld lines that gives the 
mobility to sense potential fluctuations. Thus the fluctuations are felt relative to 
the flux surface that contains the fleld line and the effective potential is ip — (f. 

The normalization used is as follows: 

t — ^phys'f^th/ Lt 

X 2^phys/ P 
y = 2/phys/p 



V = t'phys/ ^'th 



h=h 



phys" 



np 



qLt 

^ = V'phys^ 
= -Pophyst^th/^ — (27r)3/2 



where p — v^h./^, t'th = {T/f^Y^'^, ^ — qB/m and L^^ — dlnT/dx and dimen- 
sional quantities have subscript 'phys' (but only where it is necessary to distin- 
guish them with the dimensionless versions). The temperature ratio is deflned 
Te = 7: = Z'^ and Z = g^/e. 

Note that quantities are normalized in accordance with the turbulent species, 
dropping subscripts. For example, it should be understood that ITG length scales 
are normalized to pi and Lt^ in the perpendicular and parallel directions while 
ETG has normalization relative to pe and Lt^. Also it is worth making explicit 
that here — QeB/mg — —eB/rrie while = qiB/rrii — ZeB/rrii. (This 
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convention gives identical normalized equations for the linear theories of ETG 
and ITG). 



3.2.1 The primary mode 

The instabilities considered here are the ITG/ETG slab and toroidal modes. As 
described in Ref. (IGKSQip . the slab instability can be physically understood in 
terms of phase differences in the perturbed quantities that lead to an unstable 
feedback process. The slab mode has variation on the y and z directions. It is 
a three dimensional instability since it relies on a temperature gradient in the 
radial direction, mode structure and propogation in the poloidal direction and 
flows along the magnetic field. The toroidal or curvature-driven mode is much 
simpler as it is unstable in the two dimensional limit and can be understood in 
terms of equilibrium drift (v^) reinforcing density and pressure perturbations. 

Returning to local gyrokinetic theory, we write the general primary mode as 
a plane wave with variation in the y and z directions. The form of the solution 
is assumed to be 

hp = (fpo fp{vii,v±) expi{kpY + ki\z - Upt) (3.4) 
'■Pp = '-Ppo exp i{kpy + k\\z - LUpt) (3.5) 



where (ppQ is a real amplitude. Substituting this into the gyrokinetic Eqn (3.1) 
yields an expression for the kinetic response fp 



fp = JoikpV±)Fo 



kpiv'^ + vy2 - 3/2) - Up 



fcpeT(vjj + v]_/2) + f - iUp 
where the approximation a; ■ b • Vb ^ x ■ V In 5 allows one to write 



(3.6) 
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vd ■ y = erivfi + v\/2) 

The quantity is interpreted as L^/R when considering modes locahzed to 
the outboard midplane. The dispersion relation is obtained by combining Eqn 



(3.6) and Eqn(3.3) noting the flux surface average is zero for the primary mode 



solution. 



2n 
1 + r 



POO POO 

/ / MkpV±) fpV±dv±dv\i (3.7) 

Jo J ~oo 



Obviously, Eqn. 3.7 must be solved for the complex frequency ujp. One ap- 
proach is to employ a numerical search algorithm over the complex plane as a 
two dimensional space. However, since the equation does not depend in general 
on two variables, but on the single complex number ujp, an improved route is 
found in MuUer's method (see (IFraSSP ). MuUer's method is used to find complex 
roots of transcendental equations using quadratic interpolation. The convergence 
depends on the transcendental equation being analytic in the neighborhood of 
the root. The quickest method employed by the author has been to initialize 
the search algorithm to a solution from fluid theory and use MuUer's method 
to locate and polish the solution to the kinetic dispersion relation. Growth rate 
plots can then be obtained quite efficiently by making small adjustments to the 
parameters of the primary mode and using the previous solution as an initial 
guess for the new root. 

An example from a Mathematica™ calculation is given below. It is remarkable 
how closely the solution to the local dispersion relation agrees with exact (full 



geometry) eigenmode calculations. Fig 3.1 compares numerical solutions of this 



dispersion relation with results from linear runs^^^ on GYR0™1 ^j^h the 
equilibrium corresponding to the Cyclone [ncco6|i case. 
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Figure 3.1: Primary growth rate comparison: Taking rj = 3.14, = .14 
and r = 1, the local dispersion relation is solved and compared with results from 
hnear run^^sQl from GYRO.EwQl (color figure) 
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It is worthwhile to note that this hnear solution to the gyrokinetic equation 
is actually an exact nonlinear solution - linear modes are not self-interacting in 
the local approximation. To recover the fluid limit, the kinetic response given in 



Eqn. 3^ can be expanded in the slab and toroidal limits. In the slab limit we 
take er = (low-beta limit) and f ||A;|| Up <^ kp{rj~^ + t>^/2 — 3/2) (strong tem- 
perature gradient limit). This leaves gaussian integrals which can be performed 
analytically. For the toroidal branch we set k\\ = and kpexiv^^ + v'j_/2) <^ ujp <^ 
kp{rj~^ + "^^72 — 3/2) so that the temperature gradient again dominates but is 
made unstable by bad curvature. 



3.3 Introduction to secondary instability theory 

The instabilities that drive turbulence and transport in tokamaks themselves be- 
come unstable at finite amplitude to secondary instabilities. These "secondaries" 
are a key part of the nonlinear physics. The following sections in this chapter 
present a fully gyrokinetic secondary instability theory for electron temperature 
gradient (ETG) and ion temperature gradient (ITG) driven turbulence, drawing 
from the author's work ( IPluOTp . We continue in the local electrostatic limit to 
derive an integral equation for the secondary mode. The assumptions made in 
deriving the equation are designed to find "fast" secondary modes which satisfy 
7s ^ 7p and can therefore lead to mode saturation. Finite Larmor radius (FLR) 
and other kinetic effects are treated exactly capturing A;p > 1 as well as fcp ^ 1 
quasi-singular behavior. 

Aside from their potentially important role in anomalous transport in fusion 
plasmas, secondary instabilities are of theoretical interest as a mechanism for the 
transition to turbulence and as a simplification of the fully nonlinear state of 
plasma turbulence. 
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3.3.1 Motivation 



ITG turbulence is characterized by a self regulating process whereby zonal flows 
{E X B flows due to potential perturbations varying only in the radial direction) 
are strongly excited leading to the shearing of turbulent eddies which both in- 
hibits transport and reduces the linear growth drive. Internal transport barriers 
(ITB) are characterized by the presence of spontaneously created sheared flows. 
Although the role of turbulence in generating these flows is debated, they do have 
the effect of suppressing ITG turbulence and reducing ion thermal transport to 
neoclassical levels. However, within these barriers, the electron thermal trans- 
port persists at anomalous levels with a temperature gradient close to the critical 
gradient of ETG linear theory. 'ESSHi 

The presence of persistent high electron thermal transport in both experi- 
ments and simulations has been a challenge to explain. Initial mixing length 
estimates for ETG driven turbulence gave Xe ~ ^/rn^Jroi Xi which caused many 
to believe that electron thermal transport due to ETG should be negligibly small 
compared to the transport from ITG. These estimates were motivated by the no- 
tion that ETG and ITG turbulence are nearly isomorphic and might be expected 
to exhibit analogous spectrums - i.e. that typical eddy size and turnover times, 
etc., could be obtained by simple transformation from ITG values. With the 
discovery of streamer structures in ETG simulations, ''^^KQqIi j^. bg^ame clear that 
fully developed ETG turbulence is qualitatively different and has the potential 
to deliver non-negligible thermal transport. 

Mixing length theory is used to estimate transport in plasma turbulence. It 
does so by assuming a representative eddy (e.g. an ETG streamer) which stirs 
plasma in the presence of an equilibrium temperature and density gradient. There 
are questions which naturally arise if we are to interpret ETG or ITG turbulence 
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in this framework. What typical temporal and spatial scales should one expect? 
Do primary mode structures exist which are stable or quasi-stable? If so, do 
they have an important role in transport? At what amplitude do primary modes 
'saturate' and by what mechanism? In other words, what are the properties 
the "typical" eddies which cause transport? From this information, one can 
obtain estimates and the scaling properties for transport ''^2^^. We will show 
how secondary instability theory can be used to obtain these estimates and also 
gain insight into underlying physical processes. 

To get a sense of what is meant by "secondary instability" , one may consider 
the following idealized lifetime of a primary mode. It is broken into three stages 
as follows. In the first stage the mode is growing linearly. The second stage 
begins when the amplitude of this primary instability gets large (along with its 
associated gradients) and a secondary mode is destabilized with a growth rate 
substantially larger than the primary growth rate. The secondary growth rate is 
proportional to the amplitude of the primary and its time dependence is therefore 
commonly described as the exponential of an exponential. Despite the fact that 
the growth rate of this secondary is large, the perturbation itself is still small and 
therefore the primary mode growth is basically unchanged. In the third stage, the 
secondary mode has rapidly "caught up" to the primary and the dynamics become 
fully nonlinear leading to a "turning over" of both the primary and secondary 
modes. 

The secondary instability described in stage two is the subject of the remain- 
der of this chapter. The true non-linear evolution of an ITG/ETG mode is not 
as simple as the above description. Coupling with other modes is always present 
to some degree. It is therefore only in the most controlled simulations where one 
actually encounters a single linear mode growing in a quiescent background. This 
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scenario is, nonetheless, of fundamental interest. In studying secondary instabil- 
ity theory, we hope to find dominant effects which determine the features (scales, 
amplitudes, etc.) of fully developed turbulence. 

Here we consider what will be referred to as "fast" secondaries which are 
normal modes that satisfy the condition \uJs\ ^ \ujp\ where Ug and Up are the 
secondary and primary mode frequencies. In particular, the modes of interest 
also satisfy 7^ ^ 7p. Such modes are desirable to study because they are readily 
observable in simulations and are regarded as an important mechanism to bring 
about mode saturation because they have the ability to "catch up" to the pri- 
mary. Thus the primary mode functions as a "frozen" equilibrium which drives 
the secondary. The key assumption leading to the destabilization of a "fast" 
secondary is that the gradients of the primary mode (pressure, density, etc.) 
become so large that they dominate over the background gradients. Note that 
this does not violate the assumptions of gyrokinetics (which require for example 
6n <^ uq) since the spatial scales of perturbed quantities are formally an order 
smaller than the background scales allowing perturbed and equilibrium gradients 
of comparable magnitude. In the presence of such a strong primary, the effective 
total gradient becomes oriented in the direct ion ''^^^^ (on the outboard mid- 
plane, this would be in the poloidal direction for the most unstable linear modes) 
instead of the radial direction. Concurrently, the electrostatic perturbation is 
causing E x B and parallel flows which are sheared and can give rise to a Kelvin- 
Helmholtz type instability. Thus there are several ingredients present which can, 
given the appropriate primary mode phasing, create a strong secondary insta- 
bility. The assumption of large perturbative gradients also results in a simple 
version of the gyrokinetic equation where the background linear drive terms do 
not enter. (Also, geodesic acoustic modes will be absent as a result of ignoring 
background gradients.) 
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The interest in secondary instability theory for tokamak turbulence has been 
strong in recent years. Authors have explored both "fast" IBE^^ | jdko o] | cks91 |i 
"slow" ™l EESa EM CHqD (I^^i ^ |^^|) instabilities of primary mode 
structures. The break up of primary mode structures via secondary instabilities 
has been shown by many authors to be tied closely to zonal flow generation (see 
for instance Refs ( ICDOlt ICLWOO[ IRDKOOP ). leading to a predator prey- type un- 
derstanding of ITG turbulence. ''^^^^ Secondary instability theory has also been 
used successfully to estimate saturation amplitudes and explain the qualitative 
differences between ETG and ITG driven turbulence. ''^^SSli ^hus secondary insta- 
bility theory has enjoyed success in advancing our fundamental understanding of 
tokamak turbulence. 



Refs (I.TDKOni K^KSQll l(;Dn4l l(::LWnnj) investigate modes destabilized by both 



linear and nonlinear turbulent structures. We broadly refer to these modes as 
secondary instabilities. They employ a variety of model equations and ordering 
assumptions but have a common aim of investigating instabilities which arise from 
the presence of primary mode structures. The secondary instability considered 
in the present work is closely related to the one considered in Ref (IJDKOOp and 
to the secondary instability described in Ref (ICKSQip . It can be thought of as a 
fully gyrokinetic extension of those works. 

3.3.2 Content 



The remaining content of the chapter is as follows. In Sec (3.2) the equations for 
the secondary instability are derived. The primary mode is calculated using the 
local kinetic dispersion relation (which is justified for large range of perpendicular 
wavenumbers) . The amplitude of the primary is assumed to be large enough (as 
described above) that the secondary instability is driven only by the primary. 
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Therefore, the background equihbrium only enters in the secondary calculation 
through the dependence of the primary mode solution on the parameters = 
ZTe/Ti, r] = Ln/Lr and ex = Lt/R. 



In Sec (3.4) results are discussed for the case of the toroidal branch primary 
mode =0). At small primary wavenumber, the maximum ETG and ITG 
secondary growth rates are proportional to kp and kp respectively in agreement 
with Ref (IJDKOOp . The secondary instability dependence on the parameters 
1] = Ln/ Lt and ex = Lt/R is studied. The parametric dependence of the ETG 
secondary instability is qualitatively different and significantly more sensitive 
than the ITG case. 



Sec (3.5) continues the discussion of the toroidal case. Mixing length theory 
is used to explore the consequences of secondary instability theory on transport 
in tokamaks. In both ITG and ETG, the well known Dimits shift I^bbooji ^i^ove 
the critical gradient parameter {R/LT)cTit can be seen to follow from secondary 
theory. The difference between the ITG and ETG cases suggest that stiffness 
may play a weaker role in the determination of the electron temperature profile 
due to an unexpected weakening of the ETG secondary at strong temperature 
gradient (more specifically, small ey). We draw attention to recent experiments 
in electron cyclotron heating which did not find a critical gradient in the electron 
temperature . 12i2i05|i rpj^g weakening of ETG secondaries at high temperature gra- 
dient also may be a key element in the understanding of electron ITB formation. 



Sec (3.6) presents results for the secondary instability for slab primary modes 
{er = and fcy =^ 0) initially investigated in Ref flGKS91|l . The ITG and ETG 
theories are identical, in agreement with Ref (lJD02p . The mode is strong at 
small primary wavenumber kp in contrast with the toroidal case. Thus the role 
of this instability deserves a more detailed investigation especially in light of the 
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fact that the parallel wavenumber of toroidal ETG /ITG primary modes becomes 
important in the small kp limit. 



Finally, Sec (3.7) takes a brief look at secondary instability theory applied to 
the the intermediate ITG/ETG coupling regime that has been explored in recent 
numerical simulat ions . I ^'^^S H | fdo7 |i r^^ie regime p^^ < k± < is especially 
challenging because the adiabatic (aka Boltzmann) response is no longer valid 
and the kinetic treatment of ions is dominated by FLR effects. This makes 
gyrokinetic secondary theory well suited to explore the transition and interaction 
between ITG and ETG driven turbulence. 

3.3.3 Secondary instability 

The secondary instability equation is derived by considering a small perturbation 
to the primary mode 

h = hp + hs and (p = ipp + (3.8) 
with (fs/fp ^ 1 and hs/hp <^ 1. 



The solution for the primary mode is given by Eqn (3.4 - 3.7). At next order 
one obtains a linear equation for hg and As discussed, we are looking for 
"fast" secondaries and assume the primary mode dominates over the background. 
Specifically we take Vhp ^ VFq and terms involving the equilibrium do not 
appear in the equation that follows. Therefore the secondary mode draws energy 
directly from the primary mode and can serve as a saturating mechanism. The 
time dependence {ujpt) of the primary mode is taken as a constant phase which is 
set to zero without loss of generality. We thus restrict the calculation to secondary 
modes with frequency Ug satisfying \ujs\ ^ \ujp\. This includes the most important 
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modes with 7^ > 7p that are capable of "catching up" with the primary mode. 
In summary, the primary mode serves as a "frozen" fine-scale equilibrium that 
drives the secondary instability. The gyrokinetic equation for the secondary is 



(3.9) 

The third and fourth terms on the left hand side result from the Poisson bracket 
nonlinearity in the gyrokinetic equation. The first of these represents the action 
of the primary as a perturbed gradient (i.e. with Vhp in place of V-Fq)- The 
second term represents the action of the sheared E x B primary flow on the 
secondary mode. The secondary is next written as a fourier mode in z (because 



uJs\ ^ \^p\, see Sec 3.6) and x multiplied by what will loosely be referred to as 



an eigenf unction, depending on y 

ips = <f{y) expi{k^x + k\\sZ - Ust) (3.10) 
with the gyrocenter distribution function given by 

hs = h{Y,v±,v\\) expi{kxX + kusZ — Ust) (3-11) 



Combining this with quasi-neutrality Eqn (3.3) gives an integral equation for 



(p{y)= j d\n{v^,^,^')^{y-v^{sin^-sin^')) + 5(^jTG):^v{y) (3.12) 



with 



78 



sin[fcp(j;— tlx sini9)] Jo(fcpfx)+^||fc||z~'^z 



dvn v±dv±d'dd'd' 



2tt{1+t) 

^pO 



(3.13) 



The secondary stability Eqn (3.12) is solved for (p{y) and uj^ for fixed parameters 



kp, T, T], ex, and fcy^. First assume a truncated fourier series expansion (as 
used for instance in Refs (IThe92t IFreQip ). 



Ak 



n=-M 



M 



inkpy 



(3.14) 



where kf is the Floquet exponent. Eqn (3.12) then leads to the matrix equation 



c = L c + D c 



(3.15) 



with c the column vector with elements {c-m---Cm}- The operator D carries 
through the flux surface average term for ITG. The elements of D are 



Dr 



T+1 



(3.16) 



The factors 6k and 6k^^^ simply account for the fact that non-zero kf or k\^z yield 
a vanishing flux surface average of cps- The elements of L are 
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27r f 

Lmn = / '^^ll^-L'^^-L ZmnJo{kmV ±) Jo{knV ±] 



27T 







sin(?/) Jo(/i;pf_L) + f |jA;|j 



cos{y)Jo{kpV±) 



y=ypoie 



(3.17) 



with 



km = \/kl^ {kf + mkpY 



kn = ^/kl + {kf + nkpY 

■ „i ra;z-j;||fc||2 1 



ypole = sm 



JoikpVj_) . 



The matrix approach with a truncated fourier series has the advantage that 
the and integrals are done analytically to obtain Bessel functions and the 



integration over y is handled by contour integration in the final line of Eqn 3.17 
The system is solved by finding roots which correspond to functions (f pos- 
sessing convergent fourier expansions {(fM ip as M ^ oo.) 



3.4 Secondary instability for toroidal branch primary mode 

Consider the case fcy = 0. This corresponds to a primary mode driven by mag- 
netic curvature and the temperature gradient. This mode belongs to the toroidal 
branch and its secondary instability will be referred to as the Rogers-Dorland 
(R-D) secondary instability. This is the instability considered in Ref ( IJDKOOp . 



The authors of Ref (IJDKOOp identify the physical mechanism driving this sec- 



ondary as the shear in the radial ExB fiows generated by the primary. Thus the 
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R-D secondary is closely related to the Kelvin-Helmholtz instability from fluid 
theory. This statement is justified by comparing the secondary instability equa- 
tion derived from the gyro-fluid equations to the Rayleigh Stability Equation - 
the equations differ by an additive constant due to the "ion response." ''^^^^ 
(However, the difference is significant as the conventional Kelvin-Helmholtz in- 
stability is markedly weaker than the R-D ITG secondary but stronger than the 
R-D ETG secondary.) 

We find for the case of ETG that Us is weakly dependent on the Floquet 
exponent kj and that kf acts to stabilize the secondary instability. For ITG 
non-zero kf is strongly damping as it causes the flux average term to vanish. We 
therefore set kf to zero. The same appears to be true for the parallel wavenumber 
of the secondary and this parameter is also set to zero. 



With = /cp = 0, the symmetry of Eqn (3.12) ensures that complex fre- 
quencies Uz appear in sets of four. If Uz is a solution then tu*, —Uz and —u* 
are also solutions. We can therefore confine the search for solutions to the upper 
right quadrant of the complex plane. 

Growth rate curves are computed by maximizing the secondary growth rate 



over kx- Fig |3.2| shows the primary growth rate curve as well as secondary growth 
rates for the standard (Cyclone) parameters rj = 3.14, ex = -14 and r = 1. 



3.4.1 ETG 

For ETG secondaries, we find that both the large and small kp limit introduce ki- 
netic effects that cannot be treated perturbatively and make the integral equation 
fundamentally non-local. In these cases, one introduces inaccuracies in modeling 
the secondary instability with a differential equation. 

For short wavelengths, kp > 1, the exact instability equation can be efficiently 
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(a) Growth rate of primary mode (b) Growth rate of ITG secondary mode 




(c) Growth rate of ETG secondary mode (d) Radial wavenumber at maximum growth rate 

Figure 3.2: Rogers-Dorland secondary instability: Growth rates for the 
toroidal (R-D) case are plotted as a function of primary wavenumber kp. The 
secondary growth rates correspond to the fastest growing kx mode for each value 
of kp and are given relative to the amplitude of the primary mode. 
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solved with a small number of fourier modes yielding a healthy secondary growth 
rate across the domain of instability for the primary mode. 

For modes with fc^ <^ 1 the eigenvalue ujz is approximately proportional to k"^ 
and becomes very small. Since these modes are unstable for k^ ~ kp, it follows 



from Eqn (3.13) that the growth rate maximized over k^ (see Fig 3.4) is approx- 
imately of the form 7^ oc kt, and can be exceedingly small, as also reported in 



Ref (IJDKOOp . Fig 3.3 compares the gyrokinetic secondary growth rate calculated 



from Eqn (3.12) with the growth rate computed from the gyro-fluid model, Eqn 
(2) of Ref (IDJKOOp . The gyro-fiuid model has parameters fcp, k^ and r. Here we 
take r = 1 and maximize the growth rate over k^- The gyrokinetic secondary is 
solved for Cyclone parameters and maximized over kx- The agreement is reason- 
able with the gyro-fiuid growth rate differing by about 30% from the gyrokinetic 
value. As we will see, however, the gyrokinetic secondary growth rate is sensitive 
to the parameters 77 and ex and the agreement with the gyro-fiuid model will only 
hold when these parameters are carefully chosen. 

In the small kp limit, the secondary stability equation becomes quasi-singular 
with a resonant denominator. The number of fourier modes needed to accurately 
resolve the eigenfunction (p increases rapidly as kp tends to zero and the matrices 
in the numerical calculation become large. Fortunately, it is possible to expand 
the integrals in the matrix elements in the /c^ -C 1 limit which decouples the v±_ 
integral and makes the computation less intensive. Solutions obtained from this 
approximate method agree with those obtained from the exact computation. 



Fig 3.5 shows the resolution in ky needed to give properly convergent growth 



rates. The quasi-singularity occurs as the eigenvalue uo^ of Eqn 3.12 tends to zero 
at small kp. In particular, the denominator in the kernel of the integral equa- 
tion can become small in boundary layers around the nulls of the electric field 
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Figure 3.3: Secondary growth rate comparison: For the R-D secondary, 



the growth rate from gyorkinetic secondary theory, Eqn (3.12), is compared with 



values computed from the gyro-fluid model of Ref fIDJKOOp . For the gyrokinetic 
case parameters are: t = 1, r] = 3.14 and ex = 0.14. 
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(a) kx spectrum {kp = .14) 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 




0.05 0.1 0.15 0.2 0.25 0.3 



(b) kx spectrum {kp = .35) 



(c) kx spectrum {kp = .7) 

Figure 3.4: ETG spectrum: Growth rates for the R-D secondary are plotted 
as a function of radial wavenumber for three values of kp. Parameters are t — 1, 
T] = 3.14 and = 0.14 
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Figure 3.5: Resolution in ky needed for convergence of Ug to within 50%, 10% 
and 1% of it's fully converged value for R-D ETG secondary instability, (color 
figure) 

(where flow shear is maximum). Quasi-singularity was observed in Ref (IJDKOOp 
but was found to weakly affect the secondary growth rate calculation. As Fig 



3.5 demonstrates, the gyrokinetic secondary instability exhibits more sensitivity 
of the secondary growth rate to properly resolving the fine scale structure. Such 
sensitivity underscores the need for sufficient resolution in simulations where sec- 
ondaries play an important role. 

The fine scale structure can be observed in plots of the eigenfunctions in Fig 



3.6 Here we compare eigenfunctions for different values of kp. One noteworthy 
feature of these eigenfunctions is that the sharp features appear close to the region 
where the magnitude of the shear in the E x B flow is maximum (where d'^ipp/ dy^ 
is largest in magnitude), pointing to shear flow as a central driving force of the 
instability. 
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(c) kp = .7 

Figure 3.6: Rogers-Dorland secondary eigenfunctions: Re[(y9(?/)] (black 
solid) and Im[(y9(y)] (grey-solid) are plotted for the fastest growing modes. The 
primary mode ipp = ippo cos{kpy) appears red-dashed. Parameters are r = 1, 
Tj = 3.14 and = 0.14. (color figure) 
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3.4.2 ITG 



The R-D secondary instability for ITG is easily computed for all values of kp 
across the range of primary instability. The eigenvalue Uz tends to a constant 
in the limit kp —>■ and the secondary instability equation is non-singular. As a 
result, the secondary growth rate has dependence 7^ oc kp for small kp, as can be 



seen in Fig 3.2 Additionally, the value of k^ at which 7^ is a maximum {k^ 



also tends to a constant. Of course this means that for kp <^ 1, the secondaries 
have k^^rnax ^ kp so that long wavelength modes tend to be highly unstable to 
short radial wavelength perturbations. Fig |3.2K d) plots k^^max as a function of 
kp for the Cyclone parameters. In contrast to ETG, ITG secondaries have the 

is an indication that ITG turbulence should develop 
short radial wavelength structure in contrast to the "streamers" found in ETG 
turbulence. However, one cannot be sure that the fastest growing secondaries 
may are the most effective mechanism of mode saturation and it is possible that 
"tertiary" instabilities''^^^^ play a role in determining which secondary modes 
dominate the saturation process. 

3.4.3 Parametric dependence and contour plots 

The qualitative features of the R-D secondary instability described in the previous 
section hold across the range of parameters 77 and investigated in this work. 
Namely for small kp, 7^ oc k^ for ETG and 7^ oc kp for ITG and the secondary 
is unstable for all kp such that 7p > 0. Parameter scans across the rj-eT plane 



give quantitative growth rate comparisons. Fig |3.7| plots the growth rate of ETG 
and ITG secondaries across the rj-eT (t = 1) instability region {'jp > 0) for a 
fixed value of the primary wavenumber, kp = 0.21. This value of kp was chosen 
to correspond to the typical poloidal wavenumber of streamer structures in ETG 
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r] 1.^34. 

(a) ETG (b) ITG 

Figure 3.7: Rogers-Dorland secondary parameter scan: The growth rate 
Is/'^po is plotted for r = 1 and kp = 0.21. The red hne gives the approximate 
marginal stability contour and the black region outside this contour corresponds 
to 7p < 0. Also at the top of the plot the primary growth rate goes to zero at the 
critical marginal stability parameter which is roughly eT,crit ~ 0.3 (as calculated 
from the local kinetic dispersion relation). Secondary growth rates are calculated 
for the region to the right of the red contour where 7p > 0. The color scales 
linearly for ITG from black < 0.035 to white > 0.05. ETG has a more sensitive 
parametric dependence with black < 0.0003 and white > 0.003. 

simulations. The lines in these plots correspond to curves of constant 'js/fpo- The 
qualitative differences between ETG and ITG are apparent: the ITG secondary 
growth rate is generally bigger close to marginal stability of the primary mode, 
whereas the ETG secondary appears strongest at small e^. Also we note that 
the strength of the ETG secondary varies by an order of magnitude whereas the 
ITG growth has a less sensitive parametric dependence. 



From the plot of the maximizing radial wavenumber kx^max in Fig 3.8, it can 
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(a) ETG (b) ITG 

Figure 3.8: Rogers-Dorland kx^max parameter scan: Contour plot paramet- 
ric dependence of kx,max, the value of corresponding to the fastest growing 
secondary. Here we have kp = 0.21. The red line gives the approximate marginal 
stability contour and the black region outside this contour corresponds to ■jp < 0. 
Also at the top of the plot the primary growth rate goes to zero at the critical 
marginal stability parameter which is e^crit ~ 0.3 (as calculated from the local 
kinetic dispersion relation). 

be seen that the enhancement of growth rate at small for ETG is accompanied 
by a qualitative change in the radial wavenumber (k^) spectrum. For large e^, 
the instability condition k^ < kp is approximately true. But in the enhanced 



growth rate region there are unstable modes with k^ > kp. Furthermore Fig 3.8 



shows there is a small region (white) where the fastest secondary modes satisfy 
kx,ma,x > kp. This may mani 
for systems in this regime. 



kx,ma.x > kp. This may manifest as a qualitative change in the nature of turbulence 
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(a) ETG (b) ITG 

Figure 3.9: Effective aspect ratio parameter scan: Constant contours of tlie 
effective aspect ratio at saturation C,sat/^p are plotted. The primary wavenumber 
and temperature ratio are field constant: kp ~ 0.21 and r = 1. As expected, 
ETG exhibits the signature large aspect ratio "streamers," while the ITG aspect 
ratio is kept below unity aspect ratio. Both cases indicate a strong suppression 
of turbulent transport near the linear marginal stability cutoff {LT/R)crit (top of 
plots) giving a theoretical understanding of the Dimits shift. 
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3.5 Transport and mixing length theory 



Although ITG shows a much more modest variation in the secondary growth 
rate shown in Fig |3.7K b), in both the cases of ETG and ITG, it is interesting 
to consider how the variation can effect transport. To this end we consider the 
mixing length estimate of diffusivity which gives x ~ ^r/''" where is a radial 
step length and r is a nonlinear step time. 

As in Ref (IGKSQip we have in mind a typical electrostatic eddy of radial extent 



Ij, and lifetime r. Such a structure convects plasma hj E x B flow. The distance 
that the plasma can "step" is limited by the radial extent ir but also limited 
by the maximum fluid displacement determined by the strength of the E x B 
flow and the linear timescales associated with the eddy. Here we define the fluid 
displacement ^ as the time integral of the primary mode E x B velocity (which 
in our normalization is simply ve = kpipp): ^ = J kpippdt. To obtain an estimate 
for the maximum fluid displacement we first estimate the maximum amplitude of 
the primary ("saturation amplitude") as the amplitude at which 7^ = ^p.CEKQQli 



From Eqn 3.13, the secondary growth rate may be written 7^ = 7sV?po- Then 
the saturation amplitude is ipsat = Ip/ls- Finally the fluid displacement at mode 
saturation is 

^sat = ^i^p = Vsat) ~ -^^^ (3.18) 

This quantity can be interpreted within the mixing length framework as a 



radial step length. The growth rates of Figs 3.7 a) and |3.7[b) are reinterpreted 



in terms of ^sat in Fig 3.9 Constant contours of the quantity kp^sat/'^T^ = ^sat/^p 
("effective aspect ratio") are plotted in the r] — ex plane. 

Note that the color scale is logarithmic and the variation in the effective aspect 
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ratio is quite large. In both mechanism for the nonhnear shift away from 

marginal stability (Dimits shift ''^^SSSJ^ jg evident. At the top of the plot, the 
growth rate of the primary goes to zero at a critical value of er- It follows that 
the secondary growth rate can exceed the primary growth rate at small primary 
amplitude close to this marginal stability boundary. The result is a very small 
effective step length for the turbulent mixing. Therefore, transport would not 
be expected to "turn on" until the temperature gradient is shifted some amount 
above marginal stability. 

A surprising feature of the ETG case of Fig |3.9K a) is the reduction of the 
effective aspect ratio at < 0.1 where the temperature gradient gets very large. 
This is due to the strengthening of the secondary instability. Coming back to Fig 
3.9[ a), one may imagine a scenario where heating is causing a steepening of the 
temperature gradient. Keeping the density gradient and major radius constant, 
the local parameters would tend toward small er and increasing 77. Initially close 
to marginal stability, turbulent transport is effectively turned off and heating 
pushes the temperature gradient up easily. However, the effective radial mixing 
length C,sat plateaus at around ~ 0.12. After this point, ^sat now actually de- 
creases as the temperature gradient climbs. Such behavior supports experimental 
results that the electron temperature profile is not strongly stiff. ''^^^^ Also, note 
that values of less than 0.1 are observed in electron internal transport barriers 
(see for instance Ref flSGS99p or Ref (IMEJ03p l The strengthening of the ETG 
secondary at low may therefore be a key element in the formation of electron 
ITBs. 
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Figure 3.10: Cowley (slab branch) secondary: r] = 3.14. The primary growth 
rate is maximized over k\\. The secondary growth rate is maximized over both 
and 

3.6 Secondary instability for the slab branch primary mode 

For most of the spectrum of unstable wavenumbers kp, the toroidal branch is the 
dominant primary instability branch for ETG/ITG tokamak turbulence. How- 
ever, the slab branch may play a role in long wavelength turbulence'* ^™^ *' , espe- 
cially at wavelengths where the toroidal branch is stabilized. The primary mode 
with ex = but k\\ ^ belongs to the slab branch of ETG/ITG. It is character- 
ized by sheared parallel flows which are unstable to a secondary instability first 
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studied by Cowley in the fluid limit ''^SSSB^ this reason it is referred to as the 
Cowley secondary instability. 

Because the primary mode now has z-dependence, the eigenfunction (p should 
formally now have both z and y dependence. This complication can be avoided by 
noting that the 2;-dependence of is unimportant unless fthC^lnv^s/c^z ~ O{ojs)- 
The assumption that uJs ^ oJ-p then guarantees that dXv^^sjdz ^ k^^ is satisfied. 



The equation to be solved is again Eqn (3.12). The flux surface average found in 
the ITG equation is now zero and the ITG and ETC secondaries are the same. 
This is in agreement with simulations by Jenko and Dorland''^^^ that showed 
slab branch driven turbulence having the same normalized heat flux for ETC and 
ITG. 



Fig 3.10 shows results for r] = 3.14. The growth rate 7^ is maximized over 
both kx and fcys in the plot. Like the ITG R-D secondary, the Cowley secondary 
has dependence 7^ oc kp for small kp and therefore is much stronger than the 
ETC R-D secondary. The fastest growing secondaries also have finite k^^max in 

it can 



the limit kp ^ and have k\^s^rnax oc kp. By comparing the plots of Fig 3.10 



also be seen that the assumption Us ^ ujp does indeed guarantee that k\\s ^ k\\ 
and the Fourier mode dependence in z is justified. (An amplitude ippo which gives 
0;^ ^ Up also gives k\\s ^ k\\.) 



3.7 ETG secondary with kinetic ion response 

Recent simulations indicate that the adiabatic ion response (assumed so far in 
this analysis) may not capture the physics needed for resolving the non-linear 
behavior of long wavelength ETG turbulence. I ^'^^O H | fdo7 |i p'jQ turbulence simu- 
lations conducted on independently developed codes employing the adiabatic ion 
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response exhibit a divergence of the electron heat flux at high normalized shear. 
This increased heat flux is associated with the secular growth of long wavelength 
fluctuations. IM^coeji inclusion of a non-zero kinetic ion distribution function has 
been shown to give convergent electron heat transport. It is possible to examine 
this issue from the perspective of gyrokinetic secondary theory because the FLR 
effects which dominate the kinetic ions can be treated exactly. 

The mode investigated here is the R-D Secondary. The dispersion relation 
is modified to include two kinetic species and the gyrokinetic equation for the 
ions is normalized relative to the ETG scales with the exception of velocities in 
the ion kinetic response which are normalized to the ion thermal velocity. The 
secondary instabihty equation follows as before: 



(p{y) = J '^^^ '^e ~ ^±(sin?? — sin??')) — Tij Lp{^y — t>_L/i(sin'(9 — sin??')) 

(3.19) 

where 



^ sin[fcp(j/— tlx sini?)] Jo(A:pDj^)+i)||A:||2— tiJz 

'* sin[fcp(j/— t)x/isini9)] Jo(fcpiixM)+i'||fc||z— i^z 

yU = \pilPe\ = Zy^miTi/{Teme) (3.20) 

Gyrokinetic simulations capturing both electron and ion scale turbulence are 
made possible by setting the larmor radius ratio /i to be smaller than the actual 
physical value (e.g. 20 — 30 instead of 40 — 60). This technique yields a vast 
improvement in computational efficiency by essentially squeezing electron and 
ion scales together and reducing the range of temporal and spatial scales that 
need to be resolved. 
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Here we calculate (again with Cyclone parameters) the effect of including the 



full kinetic response on the secondary instability growth rate. Fig 3.11 plots three 
cases: the ETG R-D secondary with adiabatic ion (ETG-ai) response, the kinetic 
ion (ETG-ki) case with the unphysically small /i ~ 30 and the ETG-ki case with 
the physical (corresponding to deuterium) value /i ~ 60. Although the data for 
the ETG-ki and ETG-ai do not overlap, both ETG-ki cases should agree with the 
ETG-ai plot well before kpPe = 1. (The approximate continuation of the /i = 30 
ETG-ki case is drawn in dashed red and labeled "trend projection.") Starting 
from the right hand side of the plot, we see the familiar weak ETG-ai secondary 
that goes like k^. Then there is a recovery of the secondary growth rate due 
to the kinetic response of the ions. However it is clear that the value of fi has 
a dominant effect on the strength of the kinetic ion secondary (relative to the 
adiabatic ion secondary). 

The simulations of Ref (1WCF07P are invariant under scaling of /i in the sense 
that the electron thermal diffusivity measured in ion Gyro-Bohm units does not 
change significantly for several cases (/i = 20, /i = 30, fi = 40). However, 
electron thermal transport at ETG scales does increase sharply with fi when 
measured in electron Gyro-Bohm units. We interpret this to indicate that ETG 
driven turbulence becomes more intense as its coupling with ITG turbulence 
weakens and as the kinetic-ion secondary weakens (relative to ETG scales)-but 
this is accompanied by an overall downward scaling ETG transport as the ETG 
spatial scales shrink relative to those of the ITG turbulence (observe that XcBe ~ 
XGBi/f^)- So the net result is an invariance of electron thermal transport in 
absolute terms (ion Gyro-Bohms). 

We can gain insight into the role of secondary instabilities in long wavelength 
ETG turbulence by comparing the kinetic secondary instability results with Fig 
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Figure 3.11: Rogers-Dorland secondary instability with kinetic ions: Pa- 
rameters are r]i = rje = 3.14, e^i = eye = -14 and r = 1. (color figure) 
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Figure 3.12: Spectrum of electron energy transport versus wavenumber. (a) 
Large-box ITG/TEM-ETG coupled simulation; (b) small-box ETG-ki uncoupled 
(no ITG cascade). Reproduced with permission from Ref (IWCFOTp . (color figure) 



3.12[ b) taken with permission from Ref (IWCFOTp . There is a rough correspon- 



dence in the location of the trough of the secondary growth rate curve and the 
peak of the "small box" ETG turbulence spectrum. In this case it seems plau- 
sible that the secondaries play a role in determining this spectrum. (Actually, 
the quantity being plotted in Fig 3.12[ b) is thermal diffusivity per mode and 



should be compared with something like a mixing length estimate of diffusion 
using isatikp) as the radial step length-in this case the error of the location of 
the peak is still "rough," at about 50%.) However Fig 3.12K a) shows that the full 



ETG-ITG coupled simulation exhibit a monotonically decreasing transport spec- 
trum that is dominated by ITG wavelengths and does not show special behavior 
at the location of weakest secondary growth rate. This seems to suggest another 
mechanism is important for bringing about saturation of long wavelength ETG 
modes. 
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A prevalent view is that the growth of long wavelength ETG modes is mod- 
erated by ambient shear flow resulting from the ITG turbulence cascade-and 
therefore we might conclude that the kinetic-ion secondary instabilities do not 
play a critical role in the full ETG-ITG coupled scenario. 

However, one should not count out kinetic ion secondaries entirely. Indeed, 
it is useful to remember that one of the original motivating factors to giving 
attention to ETG turbulence was the persistence of electron thermal transport in 
cases where ITG turbulence was supressed'^^^^^ (i.e. internal transport barriers). 
In such a case we may expect these kinetic ETG secondary instabilities to be quite 
important and lead to a transport spectrum like the one described in Fig 3.12[ b). 



3.8 Conclusion 

Secondary instability theory is a key part of the nonlinear physics of tokamak 
turbulence. This work has shown how the local electrostatic gyrokinetic theory of 
secondary instabilities has bearing on important topics in the theory of tokamak 
turbulence. 

First we have seen that the gyrokinetic theory compares well with the gyro- 
fiuid model, yielding and kp growth rate spectra for ETG and ITG R-D sec- 
ondaries. At wavelengths significantly shorter than the electron larmor radius, 
the gyrokinetic secondary growth rate is healthy. (This eliminates the possibility 
of stable ultra-fine-scale streamers.) Quasi-singular behavior of the secondary 



stability equation Eqn (3.12) for long- wavelength primaries (small kp) calls for 
high poloidal resolution for an accurate secondary growth rate. This is relevant 
for simulations where resolving saturating mechanisms is important. 

We have seen that estimating saturation amplitudes from secondary growth 
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rates can yield estimates for radial mixing lengths and turnover times which can 
be used to estimate transport. The parametric (ey and 77) dependence of the 
radial mixing length (or "effective aspect ratio") yields a novel way to estimate 
the effect of varying local equilibrium scale lengths {Lt, Ln, R) on transport. 
These parameter scans suggest that the suppression of turbulent transport close 
to marginal stability (Dimits shift) can be understood in terms of the relation- 
ship between secondary and primary instability growth rates. It is worthwhile 
to emphasize that the parametric sensitivity of the gyrokinetic secondary is a 
somewhat unexpected result. It indicates that the gyrokinetic secondary is not 
driven exclusively hj E x B shear flow and is not as closely related to the Kelvin- 
Helmholtz instability as the gyro-fiuid secondary of Ref (IDJKOOp . 

One of the central results of the work (lPlu07p is the surprising strengthening 
of the R-D ETG secondary at small e^. First this suggests that the phenomenon 
of profile stiffness may play a weaker role for an electron temperature profile 
controlled by ETG turbulence than the an ion temperature profile controlled by 
ITG turbulence. A strong ETG secondary at small ex also immediately suggests 
the possibility that secondary instabilities play a role in the formation of electron 
ITBs. However, there are questions which remain. Small or reverse magnetic 
shear is known to be a important ingredient in the formation of electron ITBs. 
An issue not addressed in this work is how magnetic shear effects secondaries. It 
is known that the gyro-fiuid secondary is modestly destabilized by shear. 'E^EoSli 
However, since the effect of shear on gyrokinetic secondaries is not included, 
further work would be useful to determine their role in electron ITBs. 

Finally, we calculate the ETG secondary employing a kinetic ion response. 
This calculation cannot be done in the fluid limit because the ion gyro-average 
dominates the ion kinetic response. Comparison of results with coupled simula- 
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tions suggest that kinetic secondaries are an important mechanism for saturation 
of long wavelength ETG modes in cases such as ITBs where the ITG cascade is 
quenched. 
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CHAPTER 4 



The phase space cascade in two dimensional 

gyrokinetics 

The paper on which this chapter is based has been substantially updated (with 
some errors corrected). Please find the full article online: arXiv: 0904 -0243 

4.1 Introduction 

The gyrokinetic system of equations flTH68l [RF681 [HTTTl iTLSOl lPB8T| IFUSl 
IBH07P are used to describe plasma dynamics on time-scales much larger than 
the ion Larmor period. This system is essentially a kinetic theory of charged 
rings-rings formed by the fast Larmor motion of particles around the mag- 
netic field lines. These equations have been developed by the magnetic fusion 
community to describe plasma turbulence which causes the transport of heat 
and particles in fusion devices. As demonstrated in the series of recent works 
( ]SCD07t IHCD06( IHCDOSp . gyrokinetics is also appropriate for a wide range of 



astrophysical plasmas. A review of the theoretical framework of turbulence in 
magnetised plasmas is given by ( ]Kro02p . 



In this chapter, we study the simplest gyrokinetic system-the two dimen- 
sional electrostatic gyrokinetic system driven by a statistically homogeneous and 
isotropic fluctuating force. The most important scale in this problem is the Lar- 
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mor radius scale. For electrostatic fluctuations on scales much larger than the 
Larmor radius, all particles move together with the same E x B drift and a 
fluid description is correct. At turbulent scales at the Larmor radius scale and 
smaller, particles of distinct velocities have different effective E x B drifts (as 
they sample different regions of electric field during their rapid Larmor motion) 
and a kinetic description is required. At such scales, gyrokinetic turbulence ex- 
hibits strongly kinetic behaviour in phase-space due to the presence of so-called 
nonlinear phase- mixing (jDH93j). Thus fluctuations in the distribution function 
nonlinearly cascade to create fine scale structure in velocity-space in addition 
to the the two real-space dimensions-and we may think of velocity-space as an 
additional dimension to be treated on equal footing to position-space, although 
in this case having a somewhat subsidiary role to position-space as we shall see. 

In this chapter, we borrow traditional methods from fluid turbulence. This is 
made possible in part due to the fact that the gyrokinetic system bears some re- 
semblance to equations of incompressible fluid turbulence. The nonlinearity has a 
nearly divergence-free convective velocity (the E x B drift velocity due to fluctu- 
ating electrostatic fields) while dissipation occurs via a coUisional operator which 
acts via second-order derivatives (see the appendix of ISCD07P in phase-space 
and, therefore, acts on fine scales in analogy to viscosity in fluid turbulence. In 



the proper limit (see section 4.4), the electrostatic (that is, magnetic fluctuations 
are neglected) gyrokinetic system can be reduced to the Charney-Hasegawa- 
Mima (CHM) equation ( IChaTH IHMTSp . or with additional assumptions reduces 
to the the inviscid vorticity equation (the relationship between plasma dynamics 
and two dimensional fluid turbulence was demonstrated first by ITM7ip which 
describes two dimensional Euler turbulence (lKra67[ IBat69p . Thus from a fluid 
turbulence perspective, two-dimensional electrostatic gyrokinetics can be seen as 
a simple kinetic extension of fluid equations which have been extensively studied. 
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The chapter is structured as follows. We begin in section 4^ by discussing 
the gyrokinetic model, its applicability and the specific assumptions and approx- 
imations we will be using. We then discuss the dynamical coUisionless invariants 
of the gyrokinetic system in section |4.3[ It is shown that there are two basic 
coUisionless invariants, one which will correspond to the inversely cascading Eu- 
ler/CHM energy invariant and the other, which can be chosen with some freedom, 
is related to the perturbed entropy or free energy of the system. We argue that 
this invariant must cascade to fine scales in phase-space, where it is ultimately 
acted upon by the collisional operator, and incorporated into the background 
Maxwellian to increase entropy. 



In section |4.4| we explore the relationship between gyrokinetic and the CHM 
turbulence. We derive the CHM/Euler and the viscous CHM/NS equation from 
gyrokinetics to elucidate their relationship. We also discuss the two invariants 
of the CHM/Euler equation and how they relate to the gyrokinetic coUisionless 
invariants. We find that the inversely cascading gyrokinetic invariant transforms 
continuously into the inversely cascading CHM/Euler invariant ('energy'). The 
forward cascading quantity from CHM/Euler turbulence ('enstrophy') is found 
to feed into the gyrokinetic forward cascade but composes only one part of the 
forward cascading gyrokinetic invariant. 



In section 4^, we shift attention to the non-linear phase-mixing regime k ^ 1, 
where the role of phase-space has a strongly kinetic character-this regime will 
be the focus of the remainder of the chapter. This section gives a heuristic 
perspective, beginning with a description of the nonlinear phase-mixing process 
and a phenomenological analysis based upon the work by ( ISCDOTp . 

Following phenomenology, we begin a more formal analysis of the gyrokinetic 



system. Section 4.6 introduces the statistical tools and notation that will be 
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needed, followed by a list of the symmetries of the gyrokinetic system. Sym- 
metries play an important role, and in subsequent arguments we will appeal to 



the principle of restored symmetry in the fully developed state. In section |4.6.3 
we present a derivation of exact third-order statistical results following from the 
gyrokinetic equation, in the style of ( lKol41aP and (|Yag49D. 



In sections 4.7 and 4.7.8 we derive scaling laws for various spectra. To de- 
scribe these spectra, we must first provide appropriate definitions and identities, 
tailored to the needs of gyrokinetic turbulence. Scales in real-space are treated 
in the conventional way using a Fourier transform while we find an appropriate 
treatment of velocity-space scales in terms of a zeroth-order Hankel transform, 
which has a natural compatibility with the mathematical structure of gyrokinet- 
ics. Given these definitions, we can derive approximate spectral scaling laws, for 
the forward and inverse cascades, in the limit of turbulent scales much smaller 
than the Larmor radius {kp ^ 1) and velocity scales much smaller than the 
thermal velocity {pvth ^ !)• 

4.2 Equations 

Following is a brief discussion of the equations, their normalisation and appli- 
cability, and the Boltzmann response model. An important point that emerges 
from this discussion is that the model studied in this chapter will have a wide 
range of applicability, subject to specific interpretation of the fluctuating fields 
and normalisation. Readers wishing to skip these details are advised to take 
equations 4.7 and 4.8, with definitions 4.3| and 4.9 as the model equations. 



In kinetic plasma theory, each species 's' (electrons or ions) of a plasma is 
described by the distribution function which is a function of space, time and 
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velocity. As is customary with so-called delta-/ gyrokinetics, the full distribution 
function is split into an equilibrium part and fluctuating parts such that the total 
distribution function is expressed 



/. = Fos - ipFo, + hs (4.1) 

where 9? is the electrostatic potential (the normalisation is explained below). 
The equilibrium distribution Fqs is a Maxwellian distribution in velocity-space, 
Fqs — (27r)~^/^e~^-L^^~''ii^^. The second term the so-called Boltzmann part of the 
distribution. And the final term is the gyro-centre distribution function, hs-, which 
depends on the perpendicular velocity (magnitude of the velocity perpendicular to 
the mean magnetic field), the parallel velocity (velocity parallel to the equilibrium 
magnetic field) and the gyro-centre position R which is the spatial position r 
minus the Larmor radius vector p (defined in normalised units below). It can 
be shown by rigourous asymptotic expansion that the gyro-centre distribution 
function hg satisfies the nonlinear gyrokinetic equation. In the normalisation we 
choose, the gyrokinetic equation is the same for each species. Thus we drop the 
species label. The electrostatic gyrokinetic equation, without variation in the 
background magnetic field or temperature, is 



|.„|.v..V.-(C„). = ^.„.^ (4.) 

where the gradient is in the gyro-centre coordinate, i.e. V = Vr, and Fq — 
(27r)~^/^e~''-L''^~^li'^^. We have also introduced a general kinetic forcing term 
which will be left unspecified at this stage. The normalisation is summarised in 
the following table, with physical quantities having subscript 'p': 
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t = tpVth/Ln X = Xp/p y = yp/p z = Zp/L„ 
v = Vp/vtt <^ = <^pfe h = hp"-^ Fo = Fopvl/no 



where we have the following definitions: The equilibrium density and temperature 
of the species of interest are no and Tq; the mass and charge are m and g; the 
thermal velocity is Vth = \/T\)J^s', the Larmor radius is p = fth/^c where the 
Larmor frequency is Vt^ = qB/m. The equilibrium density scale length is taken 
to be Ln = \dlnnQ/dx\^^ . The gyro-average {.)-^ is an average over the gyro- 
angle with gyro-centre position R held fixed. The real space and gyro-centre 
coordinates are related by the transformation r = R + p where the Larmor 
radius vector is p{'d) = z x v = v±{ycos'd — xsin-i?). In this normalisation, the 
gyro-average acting on an arbitrary function of position A{r) is defined 

(A(r))^ = (27r)-i I d^R + (4.3) 



Lastly, the gyro-averaged E x B velocity in equation 4.2 is 



^rE = zxV (ip)-^ (4.4) 

The full plasma dynamics are described by the evolution of the gyrokinetic 
equations for all species, with the additional constraint that charge neutrality 
must be maintained (the quasi-neutrality constraint). However, due to numerical 
and analytic difficulty, one does not typically solve the gyrokinetic equations for 
multiple species simultaneously. It is common to solve a single gyrokinetic equa- 
tion, and exploit the separation in scales between the ion and electron Larmor 
radii to obtain a Boltzmann (or 'adiabatic') approximation for the other species. 
Thus, conceptually, one can consider two separate regimes, one corresponding 
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to non-trivial ion-dynamics at the ion Larmor radius scale and ion thermal ve- 
locity scale ('ion-scale' turbulence) and the other corresponding to dynamics of 
the electrons at the electron Larmor radius scale and electron thermal velocity 
('electron-scale' turbulence). (These scales are separated in magnitude by roughly 
the square root of the ion-electron mass ratio.) 

For ion-scale turbulence, the electron gyrokinetic equation is dominated by 
the parallel streaming term which gives the simple equation v\\dhf./dz = 0. Thus 
the electron gyro-centre distribution function must be a constant in the magnetic 
surface (flux surface) defined by the equilibrium. We will make the assumption 
that this constant is zero and thus the perturbed electron distribution function 
is given only by the electron Boltzmann respons^ (normalised to ion-scales). 

Sfe ^ - V^^o (4.5) 
qTe 

As an alternative to the Boltzmann response, one may consider the 'no re- 
sponse' model, 6fe ~ 0. This may be more appropriate in the case of exact two 
dimensional turbulence because a small but non-vanishing parallel wavenumber 
is required to physically establish the Boltzmann electron response. For electron- 
scale turbulence, the ion gyrokinetic equation, expanded in the root-mass-ratio 
becomes simply dthi = 0. Thus hi is a constant which is set to zero. And the 
appropriate ion response model is the Boltzmann ion response (normalised to 
electron-scales): 

^Note that a modified Boltzmann response (see equation 2.69 of IDorOS]) would be more 
correct than this pure Boltzmann response. The modified response is physically more accurate 
as it accounts for the fact that the parallel mobility of the electron species only allows them 
to establish a Boltzmann distribution relative to the flux surface to which they are confined. 
However, the additional term in the modified response (containing a flux-surface-average of the 
potential) formally breaks the isotropy assumption but does not otherwise substantially affect 
the analysis that follows. Thus we will omit it for simplicity. 
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Sf.^-^^^Fo (4.6) 

We will now give the gyrokinetic system, normalised to the scales of the tur- 
bulent species. The subsequent analysis of this chapter will proceed without 
reference to a specific species or response model. As a final point before intro- 
ducing the normalised system, we note that although it may be useful to think 
in terms of ion-scale turbulence for comparison with the Hasegawa-Mima tur- 
bulence, the results may be even more applicable to electron-scale turbulence as 
the ion Boltzmann response becomes more exact at high wavenumbers whereas 
the Boltzmann-electron approximation will be broken for high wavenumbers - 
that is, the electrons become kinetic as the scale approaches the electron Larmor 
radius {kp^ ~ 1). 

We take equation 4.2 and neglect the parallel ion inertia term viif (see 



4.8.3 



for a discussion of this). The fy -dependence is no longer important (except for its 
appearance in the collisional operator). Thus it will be concealed when possible 
in the remainder of this chapter and it should be assumed that wherever velocity 
integration is present, integration over v\\ is impliedj^ Henceforth, we will notate 
the perpendicular velocity v± as simply v. We express the gyrokinetic equation 
compactly in terms of the gyro-centre dependent quantity g = h — {(p)^Fq. Thus 
the two dimensional gyrokinetic equation we will be using in this chapter is 



dt 



+ YE-Vg={C[h])^ + J^ (4.7) 



The gyrokinetic system is closed with the quasi-neutrality constraint, which states 
that the ion charge and electron charge, which are expressed as integrals over the 

^When written explicitly, the full collisional operator is the gyro-averaged linearised Lan- 
dau operator (see the appendix of 1SCD07]| . an integro-differential operator with non-trivial 
dependence on tiii-space. 
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distribution function, must sum to zero 



27r J vdv {g}^ = aip - Toip 



(4.8) 



where we define the angle average {.)^ as an average over gyro- angle with real- 
space coordinate r held fixed: 



{m))r = (2^) 



diM(r - p(?9)) 



(4.9) 



The operator Fq is defined Toip = f vdv e ^^/^ {{v)B)r- This operator is multi- 
plicative in Fourier-space-it is simply fo(A;) = /o(A;^)e~'^^, where Iq is the zeroth- 



order modified Bessel function. The constant a in equation 4^ varies depending 
on the choice of scales and response model 



qTe 
qTe 



ion-scale, Boltzmann- electron response 



= < (1 + ^) electron-scale, Boltzmann-ion response (4-10) 

1 no-response model 

(Note that the constant temperature ratios are determined by the fact that we 
normalise the Boltzmann response to the scales of the turbulent species.) Quasi- 
neutrality in Fourier space is expressed 



(p{k) = (3{k) I vdvJo{kv)g(k,v) 



(4.11) 



with 
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a - to{k) 



(4.12) 
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We will take the system defined by equations 4.7 and 4.8 (or equivalently equation 
4.11) to be our gyrokinetic model. We propose that this system may be thought 
of clS db simple paradigm for kinetic turbulence. 



4.3 Collisionless invariants 

In two dimensional gyrokinetics, g'^ is collisionlessly conserved for each value of 



V individually. Multiplying the gyrokinetic equation 4/7 by g and averaging over 
gyro-centre position R we have 

y^2l = y + ^ (4.13) 

where the injection rate of kinetic free energy is defined 



-^9:f (4.14) 

Multiplying this equation by F^^ and integrating over velocity-space yields 
the global balance for what we will call 'generalised free energy': 



dt J V J \ Fq 
where the generalised free energy is defined 



27T vdvl " ' ' ) +£ (4.15) 



W, = 2nJ vdv J (4.16a) 
= 27r fvdv f^^-Q^ (4.16b) 



V 2F 







and the injection rate of free energy is defined 
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Of course, multiplied by any function of v is conserved in the absence of 
collisions. (In fact any power of g multiplied by any power of v is conserved.) 



Beginning in section 4.6 , when we develop the statistical theory of gyrokinetics, we 
will re-define g to be multiplied by an arbitrary function of v so that the following 
analysis will be more general. Here, we take the preceding definition which relates 
in a straightforward way to the three dimensional invariant of gyrokinetics as we 
discuss below. 

There is a second invariant in two-dimensional gyrokinetics, which is a real- 



space invariant. It is found by multiplying the gyrokinetic equation 4.7 by (v5)j^ 
and integrating over real-space coordinate r and velocity v. The resulting equa- 
tion is 



dE _1 f (Pv 
~dt ~ 2 



J^^^^J vdv mh])^)^ + (4.18) 



where the invariant E is defined 



2JV 

and the injection rate of E is defined 



\ [V-</^ro</^] (4.19) 



1 f (Pr 

27r / vdv (J^r (4.20) 



2JV 

In this chapter we will take Wg and E to be the fundamental invariants. Ob- 
viously there is freedom to this choice. Another choice is the three dimensional 



invariant in gyrokinetics. That is, the gyrokinetic equation 4.2 retaining the term 
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v\\dh/dz, has just one invariant, just as three dimensional fluid turbulence has 
one invariant (energy). In the electrostatic limit, the single gyrokinetic invariant 
(in physical units) is the free energy, W = TqSS, where the perturbed entropy 
is 6S = - T.sl d^^{6fsf/2Fos fsee lKHQl IS()H96b|) . (In electromagnetic gy- 
rokinetics, (IHCDOGIl have derived the corresponding invariant which they refer 
to as 'generalised energy.') We can easily demonstrate the conservation of (nor- 
malised) free energy W, for the two dimensional electrostatic case by multiplying 



the gyrokinetic equation 4.7 by h/Fo and integrating over velocity and position 
r. The quantity is expressed in terms of h for comparison with (1HCD06P : 



2n I vdv^-^ 



(4.21) 



Alternatively, free energy can be expressed as a linear combination of Wg and E: 



W=W„ + E 



(4.22) 



where it may be useful to see that Wg written in terms of h is expressed 



1 f (Pr 



27r / vdv'^-^ + i^TQ^ - 2a(p'^ 



As we will see, E becomes subdominant to W„ at high-fc and we will have W ~ Wg 



in this limit. 



4.4 Charney— Hasegawa— Mima/Euler turbulence 

In this section we will derive a single equation which corresponds to the Euler 
equation or the inviscid CHM equation depending on the choice of a parameter, A. 
The system is still kinetic in that there is a passive cascade of a kinetic invariant. 
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We will derive the invariants of this system and provide a complete description 
of the relationship of these invariants to the gyrokinetic invariants derived in the 
previous section. 

At the end of the section, we will briefly show how raising the ordering of the 
collision operator can produce the viscous CHM and NS equations. Thus we verify 
that the kinetic description is needed in the regime of low-collisionality, while in 
the high-collisionality limit, the fluid equations are a complete description. 

We begin by repeating the gyrokinetic system: 



g+v^-V(7=(C[/i])j, + ^ (4.23a) 
2ti I vdv {g)^ = aip- Voip (4.23b) 



To arrive at the CHM/Euler equation from the system 4.23a, one takes the small 
ion Larmor radius limit, /c <^ 1, as well as the small temperature ratio limit such 
that (1 — a) ~ (9(/c^)|^ In this limit, the gyro-average on the E x B velocity acts 
as unity. Also, the gyro-average on the collision operator can be neglected since 
its corrections are higher order. We must keep Fq to second order in k (the zeroth 
order part gives no contribution upon substitution into the nonlinearity) , so we 
take Fq ~ 1 + in the quasi-neutrality constraint. The gyrokinetic system in 
this limit is given by 



^ + V£;o- V(7 = CM + ^ (4.24a) 
2-n I vdv (g)^ = AV - VV (4.24b) 



■^Thc limit of small temperature ratio is somewhat artificial but leads to an elegantly simple 
fluid equation, the CHM equation. We note that a more detailed system of fluid equations 
may be worked out for a temperature ratio of order unity but, for our purposes, will not be 
substantially different that the familiar CHM equation. 
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where v^o = z x Vty? is the E x B velocity (without gyro-average) and = a — 1. 
The scale given by the inverse of this constant, A~^, corresponds to the Rossby 
deformation radius in the atmosphere or the sound Larmor radius, ps, in a plasma. 
For two dimensional Euler turbulence the parameter is zero, A = 0, corresponding 



to the no-response model given in equation 4.10 Now, by integrating equation 



4.24a over velocity-space (noting that the integral of the collision operator is zero 

we 



by particle conservation) and substituting quasi-neutrality, equation 4.24b 
obtain a single equation for the electrostatic potential. This equation can be 
recognised as the inviscid and CHM equation: 

dt{\' - V^)ip + v^o ■ V(-VV) = ^CHM (4.25) 

where ^chm = 2vr J vdv {T)^. Again, the two dimensional Euler equation follows 
from this equation by taking the no-response model, A = 0. Note that there 
is traditionally a term dy(p due to the background density gradient. Although 
we have not included this term explicitly in our gyrokinetic system, it may be 
included in the source term Tc-am if needed. 

It is well known that there are two invariants of the CHM/Euler equation. 
These are referred to as energy and enstrophy, although their physical interpreta- 
tion depends upon the specific scale of interest. Energy and enstrophy, are found 



by multiplying equation 4.25 by ip and V (p respectively and integrating over the 



system volume. The invariants are 



1 f£r 

-C^CHM - 2 



j ^(AV^ + lVv^r) (4.26a) 
^CHM = \j ^(A^|Vy.p + (W)') (4.26b) 
It should be clear that the two dimensional Euler invariants are recovered by 
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neglecting the terms proportional to A^. We now show that there is an additional 
invariant associated with the (equivalent) gyrokinetic description given by equa- 
tion |4]24 This will complete our discussion of the relationship between the CHM 
(and two dimensional Euler) cascade and the gyrokinetic cascade. To derive the 



new invariant we must solve for g in equation |4.24a[ Without loss of generality, 
we introduce the following ansatz for g: 



g = F^{\^-V^)^+~g 



(4.27) 



It is clear that quasi-neutrality, equation 4.24b implies 



vdv {g)^ = 



Furthermore, by multiplying the CHM equation 4.25 by Fq and subtracting it 
from equation 4.24a| we obtain a separate equation for g: 



dt 



+ v,,o- V^= (C[^])i, + (^-Fo^c 



(4.28) 



This implies another collisionless invariant (in addition to the Wg invariant which 
must obviously still be conserved): 



1 fd^r 



(4.29) 



'^J V J Fo 

Expanding the generalised free energy invariant Wg in terms of the solution for 
g we can express it in terms of the three gyrokinetic-CHM invariants: 



Wg = Wg + X^EcKM + ZcHM 



(4.30) 
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each of which are conserved individually in the CHM/Euler limit, as we have 
shown. 

To complete the picture of how the gyrokinetic cascade continues in the 
CHM/Euler limit, we note that the the gyrokinetic invariant E reduces to the 
energy invariant with the CHM approximation Fq ~ 1 + V^. 



CHM 



(4.3i; 



As a last note, notice that in the CHM/Euler limit, the g distribution function. 



governed by equation 4.28, takes on the role of a hidden passive scalar, being 



advected by the E x B flow. 



4.4.1 Collisional limit, Navier— Stokes and viscous Charney— Hasegawa- 
Mima equations 

Viscosity does not appear in the preceding fluid limit of the gyrokinetic sys- 
tem. This is to be expected, as the collisional dissipation is ordered small in 
gyrokinetics, i.e. the collisionality is sufficiently small that a kinetic description 
of the plasma is required. To recover the viscous CHM and the two-dimensional 
Navier-Stokes equation, we must raise the ordering of the collisional operator. 
Specifically we will assume that C[g] ~ 0{k~^dg/dt). We apply the same ansatz 
for the solution of g given in equation |4.27 This time, the dominant-order equa- 
tion for g is simply 



C[g] = (4.32) 

Note that the collision operator acting on the other part of g is zero since it has 
been defined to be proportional to a Maxwellian. As the collision operator has 
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been linearised, the meaning of this equation is that colUsions between g and 
the background MaxweUian are neghgible. Consequently, g must be a perturbed 
Maxwellian, Fo{no + An, Tq + AT) — Fo(no, To). Because quasi-neutrality imphes 
J vdvg = 0, there is no density perturbation so we have g — — 3/2)FoT, 

where T = AT/Tq and T satisfies a passive scalar equation. 

At next order we obtain the CHM/Euler equation by integrating the remnants 
of the gyrokinetic equation over velocity. The integral of the collision term is 
somewhat involved. Note that we must retain all gyro-average and angle- averages 
with the coUisional term, as the dominant terms will all be zero and the final result 
is due only to the so-called finite-Larmor-radius 0{k'^) corrections. The result is 

2n J vdv{{C[{^)^Fo])^)^ = -z/VV (4.33a) 
f 

iy = - vdvd'&j{l + cos'^'&)C[v'^{l/2 + cos^-3)Fo\ (4.33b) 
Thus we arrive at the viscid CHM equation 



dtiX'' - V^)<p + veo ■ V(-VV) = -i^VV + ^CHM (4.34) 

Again, we note that the two-dimensional NS equation is recovered in the no- 
response limit, by setting A = 0. 

As we have reviewed in this section, the small-A; limit of gyrokinetic turbulence 
imparts a fiuid character to the plasma dynamics. In particular the phenomenon 
of nonlinear-phase mixing is crucially absent in this limit. In this chapter we are 
especially interested in the case where the kinetic behaviour is nonlinear in char- 
acter, which is corresponds to k ~ 0{1) and larger. In the large- A; hmit, we will 
show how statistical self-similarity and scaling relations can be derived for the 
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gyrokinetic turbulence cascade. In the following section we sketch phenomeno- 
logical arguments corresponding to the low- A; CHM limit, and then proceed to 
the phenomenology of the large-fc limit and build on the previous sections to 
give a general sketch of the full range of the stationary spectra, from CHM to 
gyrokinetics. 



4.5 Phenomenology 

We introduce the notation by sketching the phenomenology for stationary, driven 
CHM turbulence. We define to be the electrostatic potential at scale spatial 
scale I. We take this to be a fluctuating quantity-i.e. it is a random variable with 
negligible mean. We notate the averaged quantity as which may be interpreted 
as the root-mean-square potential at scale t or the square root of the structure 
function. 

We assume injection at a scale /c,, giving rise to two inertial ranges, the inverse 
cascade inertial range corresponding to k <ti ki and the forward cascade inertial 



range k':^ ki. From equation 4.25 we have 



We multiply this equation by C,~'^^pi and assume constancy of enstrophy flux in 
the forward cascade range-i.e. we take the left hand side equal to a constant, 
the enstrophy injection rate, Sz ~ T^\{^^'^\^ — £^'^)Lpj. Thus ^~^<^f ~ £z so 
that i^e ~ ^z'^^^ s-iid the yj^-spectrum is E^p ~ k~^. Likewise, multiplying 



4.35 



by ifi and assuming constancy of energy flux in the inverse cascade range gives 



~ e'J^i^/' or ~ A;-iV3. 



The cascade of Wg is passive (see for example ILHOGD and is directed forward. 
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In the forward cascade range we have e ~ £ '^gj'^e- Therefore we have 

££1/2^-1/2 _ ^^1/2^-1/6 _ 

4.5.1 The nonlinear phase-mixing range, A; ^ 1 



fISCDOTp have produced a phenomenological description of the large-fc forward 



cascade in two-dimensional gyrokinetics (see also ISCDOSp . The scaling results 



obtained by fISCDOTp are in agreement with the results of this chapter. We 



reproduce their arguments here, with an expanded notation, and also give the 
inverse cascade phenomenology. The phenomenology will motivate the following 
material in this chapter and highlight the salient physical processes. 

In gyrokinetics, ions and electrons are acted upon by fluctuating fields which 
are averaged along their gyro-orbits. Two particles sharing the same gyro-centre 
but having different velocity perpendicular to the equilibrium magnetic field are 
subject to a different averaged electrostatic potential and thus experience a dis- 



tinct E X B drift. This is illustrated in figure [41] One can see from this cartoon 
that if the E x B flow has a correlation length i, then the particle motion must 
be correlated in velocity-space such that ~ i. 

Now let Qi be the the distribution function g at scale i. Again, this quantity 
is assumed to be random and the scaling we will find will be for the averaged 
quantity gi. The central assumption for our treatment of velocity-space depen- 
dence is a two-scale dependence. To be explicit here, we may assume that g£ may 
be separated into its fine and smooth-scale dependence as follows 

ge = giF{v)Mv) (4.36) 
where F{v) is a smoothly varying function of velocity which decays strongly 
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Figure 4.1: Correlation in velocity-space 
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at large v so that we may say j v"'F"^{v)dv ~ 0{1) for any m and n. We assume 
fe{v) is a random variable which varies on the scale ~ i, as argued above, 
and has negligible mean and a variance which is of order unity. To see that ge 
functions as the amphtude of gi we may calculate the velocity integral of gj, 
exploiting the scale separation in velocity-space: 



/« OO PV + 1 

vdvgl = f, / vdvF\v)f! ^ -g^Y.'^F^v) / Jnv')dv' 
J 

~ gj I vdvF^{v) ~ (4.37) 



where we have used that JJ^"^ fj{v')dv' ~ 1. Now, from quasi-neutrality, 4.8 
we can determine the scaling of ipi in terms of gi. Equation 4.8| gives (pi ^ 
J vdvJo{v / tjgt. The right hand side is calculated as follows 



vdvJo{v/i)ge ~ i^^^ j v^^^dv cos{v/e - 7r/4)^^ 
~ f/^J^'^'/'Fiv) / Mv')dv' 



~ ig£ j v^/^dvF{v) 

~ igi (4.38) 

where we have used the large-argument approximation of the Bessel function Jq 
on the first line. Also we argue that the / integral accumulates like a random 
walk to write f^{v')dv' the result is 

~ tgi (4.39) 



From the gyrokinetic equation 4.7 we have the nonlinear decorrelation time 
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r. 



NL 



I ~ e-^ipeJo{v/£) ~ r^/\^^/V^ cos{v/£ - tt/A) (4.40) 



where we have again used the large argument approximation for Jq. Averaging 
this expression, and then multiplying by gj we get 

'-'^'v~'/'^,g', (4.41) 



Dividing by the background distribution function Fq, integrating over velocity 
and assuming a constant free energy flux gives 



r^^^ipif, ~ / vdv-^ ~ e (4.42) 



Using the result 4.39, we have ^^'^ ~ e so that we can write ~ e^/^f^/^ 



and ~ e^/^p/^. The resulting spectra are Wg{k) oc k and E^{k) oc k ^^/^^ 
where E^p is the ip^ spectrum. 

The inverse cascade range is analysed in the same manner. We integrate the 
(phenomenological) gyrokinetic equation over velocity to get 



~ £ / f Jq (t>/£)(yf£ ~ £ ^(^i / dv cos^ {v/i — tt/ 4) gi 
Tnl J J 

~ r^^^ifiigi ~ r^/2(^| (4.43) 

where we have again assumed a random walk accumulation to write J dv cos^ {v / 1— 
7r/A)g{: ~ i^^'^gi- And assuming constancy of flux of the invariant E in the in- 
verse cascade range, we have ^pI/tnl = Se ^ Thus ~ e'l^i^l'^ and 
gi ~ el^t'^l'^ and the spectra are E{k) ~ E^{k) ~ k^"^ and Wg{k) ~ k^ . 
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4.5.2 The cascade through multiple scales 



Using the phenomenological scalings from this section, we may give a sketch of the 
fully developed cascade from CHM/Euler limit to the high-fc gyrokinetic limit. 
As there are multiple special scales and invariant quantities, there are several 
possible scenarios. For simplicity we consider just two scenarios given by figures 



4.2 and 4.3 Equation 4.30 states that the gyrokinetic invariant Wg is composed 



of three separately conserved quantities in the CHM limit: Wg = Wg + X^E + Z. 
We use this relation to trace out the fate of Wg in the long wavelength limit. 



4.5.2.1 Injection at large scales 



In figure 4.2 injection is at large scales, in the so-called potential limit of the 
CHM equation, where fcj -C A and the invariant E reduces to E ^ X^E^. In the 
stationary state, the injection of E fuels the inverse cascade and the injection of 
Z and Wg feed into forward cascades which ultimately drive the cascade of Wg 
at the Larmor radius scale ~ 1. In the so-called "kinetic limit" of the CHM 
equation, A ^ A; <^ 1, the contribution of E to Wg is negligible and Wg = Wg + Z. 
In the potential limit, k <^ X, the steep spectrum of E is the dominate part of 
Wg, i.e. Wg ~ X^E. This leads to the unexpected result that there must be a 
continual accumulation of the forward cascaded quantity Wg at large scales, in 
the fully developed state. 

4.5.2.2 Injection at small scales 



In figure |4.3| injection is at small scales, in the nonlinear phase-mixing regime, 
^ 1. The free energy Wg cascades forward to the dissipation scale whereas 
the second invariant E cascades inversely. Again, the identity of Wg changes in 
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Figure 4.2: Scenario 1: Injection at large scales 
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Figure 4.3: Scenario 2: Injection at small scales 

the long wavelength limits. In the so called kinetic limit of the CHM equation, 
A <^ <^ 1, we have Wg ~ Z since there is no source of Wg in this scenario. In 
the potential limit A; -C A, we have Wg ~ X^E and the Wg takes on the steeper 
spectrum of E. 



4.6 Statistical theory 



us- 



In this section we will analyse the gyrokinetic system, equations 4.7 and 4.8 
ing familiar statistical methods from fluid turbulence (we adopt somewhat the 
methods, terminology and notation of IFriQSp . First we list the symmetries of 
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gyrokinetic system and then formulate the statistical assumptions which follow. 
Following this we derive the two exact third-order statistical relations which de- 
termine self-similarity in the inertial ranges for the forward and inverse cascades. 

4.6.1 Symmetries of the gyrokinetic system 

Statistical arguments in turbulence theory are anchored to the underlying sym- 
metries of the dynamical equations of interest. Although specific boundary condi- 
tions, initial conditions and forcing mechanisms break these symmetries at large 
scale, a guiding principle of fluid turbulence is that at sufficiently small scales, far 
from boundaries, a fully developed turbulence will emerge where the symmetries 
of the dynamical equations are restored in a statistical sense. We will appeal to 
this principle in the arguments of this chapter. 
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The following table is a summary of symmetries of the gyrokinetic system: 



Position-translation: 



Time-translation: 



Constant potential shift: 



Scaling transformations: 



Rotations: 



Parity: 



(r) - (r + a) 

(4.44a) 

(t) - (t + r) 

(4.44b) 

(g, </^) ->((? + (a - l)i^Fo, ip + ij) 

(4.44c) 

(p, r, V, t) {-f~^g, /i^7~V, Air, fiv, -ft) 

(4.44d) 

(r) ^ (Ar); A G S0(2) 

(4.44e) 

{g, ^, r) ^ (-^, -r) 

(4.44f) 



Translational invariance in time and position-space are obviously the basis for 
the assumptions of homogeneity and stationarity respectively. The invariance of 
gyrokinetics to shift by an arbitrary constant potential is somewhat trivial since a 
constant potential produces no electric field and thus has no effect on dynamics. 



(Note that for the transformation 4.44c, the shift in g is Maxwellian so as to 



vanish under the collision operator.) The scaling symmetries are described by 



the transformation 4.44d and parameterised by two scaling factors, /x and 7. 
This symmetry will be crucial in determining the scaling in the velocity-space 



increment. In particular, for 7 = 1, the scaling invariance 4.44d implies a dual- 



scaling in velocity and position-space as described in section 4.7.1 Note that the 
scaling invariance applies to the coUisionless limit. 
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4.6.2 The ensemble average 



We will denote the ensemble average with an over-bar so that the correlation 
function of g between two points in phase space, (r, v) and (r', v'), is defined 



Homogeneity in position-space (or equivalently, in gyro-centre space) follows from 
the translational invariancc of the gyrokinetic equation and implies that G is 
a function only of the increment £ = R' — R. In contrast, velocity-space is 
fundamentally inhomogeneous as there is obviously no translational invariance 
in i>-space. (Velocity dependence in the Maxwellian and the gyro-average each 
break translational invariance.) However, we expect a local homogeneity to apply 
to velocity-space for small enough velocity increment iy = v' — v . In other words, 
the two-point statistical quantities such as structure and correlation functions 
should have a two-scale dependence in velocity-space, varying strongly with the 
increment but weakly (smoothly) with the centred velocity v — {v' + v)/2. 

The remainder of this section will be devoted to deriving exact statistical 
relations in analogy to Kolmogorov's 'four-fifths' law (or the 'three-halves' law 
for two dimensional NS turbulence). As there are two conserved quantities, there 
will be two third-order statistical relations. We begin with the relation implied 
by the conservation of Wg and finish with that for E. 

4.6.3 third-order Kolmogorov relations 

Let's redefine g to be normalised by an arbitrary mean-scale function of velocity- 
space, K,{v). This will allow us to develop a theory for arbitrary velocity-integrated 
moments of the invariant g^. (We will fix k at a specific point in the analysis. 




(4.45) 
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later. This is done for the sake of clarity in presentation and we will demonstrate 
that the results are valid for general k) The definition of g we now adopt is 



g = {h- Fo{<^)-^)/k{v) 



(4.46) 



The generalised free energy, originally defined by equation |4.16a[ is also redefined. 
Using the ensemble average instead of explicit spatial averaging, we define 



Wg = 2n 



vdv 



9' 



(4.47) 



4.16a 



and we 



1/2 

which for k = Fq corresponds to the definition given by equation 
will see that this instance of Wg corresponds to the free energy W in the large-/c 
limit. 



From the gyrokinetic equation we derive the budget equation at scales i and 



ijr< 1 

^ - 2 V£ ■ S3 = {C'[h'])^, /k' + g' {C[h])^ /k + g:F'/K' + g':FlK (4.48) 



where we have defined the third order structure function ^3 = SwE^g"^ with 5g = 
g' — g and ^v^; = v^ — v^;. Also, we have used homogeneity and incompressibility 



of \e to manipulate the nonlinear flux to obtain —g'^E ■ 9 + 9^'e ' ^'9' = 

V, ■ 53/2 

Now if we take £ = 0, = 0, the nonlinear term is zero by homogeneity. The 



result (equivalent to equation 4.13) is the 'global budget' equation: 



dg^ 
dt 



2g {C[h])^/K + 2gJ^/K 



(4.49) 
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assuming stationarity, the left hand side of the equation is zero and we have the 
steady state balance between forcing and dissipation 



gJ'/K = -g{C[h])^/K = e{v) (4.50) 

where we have defined the kinetic free energy injection rate s{v). To derive the 
'inertial range' law for we now assume that the forcing is restricted to a narrow 
range around a injection scale which we denote li. For simplicity we also assume 
that the only special velocity scale of the problem is the thermal velocity v = 1 
and we therefore assume that the forcing acts on this scale in velocity-space. For 
i £i and <ti 1, the forcing term in the budget equation |4.48| reduces to 



gj^/n' + g'J^/K ^ 2gJ^/K = 2e (4.51) 

Now we take the limit of small coUisionality. That is, since the collision opera- 
tor is proportional to a normalised collision frequency (which is the collision 
frequency of the turbulent species divided by the diamagnetic frequency Vth/ Ln) 
we are free to take this parameter as small as we like so that for fixed i and 



the collisional term is negligible in the budget equation 4.48 Thus by assuming 



small coUisionality, stationarity £^ <^ 1 and i <ti ii we have from equation 4.48 



V, ■ 53 = -H^) (4.52) 

where we have written the injection rate e as an explicit function of v. Using 
isotropy, it follows that 

53 = -2l(e£ + const.) (4.53) 
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Note that the constant of integration is included here, which in general depends 
on V and In fluid turbulence, this constant is zero since the structure function 
in that case depends only on £, so must be zero when ^ = 0. However, with 
non-zero C.^ we must allow for possibility that tends to a constant as £ — > 0. 
Indeed, we have not yet determined the £t,-dependence of Sz (i.e. e does not 
depend on l.^) and we will use the additional term to account for this dependence 
in section 14.7. 11 



4.6.4 The E third-order Kolmogorov relation 

To derive the scale-by-scale budget equation for the invariant E we proceed anal- 



ogously to the analysis leading to equation |4.53[ However in this case, we are 
dealing with functions of real-space variable r and its increment = r' — r. 
Also, velocity will appear as a dummy variable of integration, so it will simplify 
our notation to refer to a single variable v in what follows-i.e. R = r — p(f), 
R' = r' — p{v) and g' = g(R', v), etc. , from which it follows that £ = £r and we 
will refer to a single increment in position-space, £. As our interest has shifted 
away from velocity dependence, let's also take the simplification k = 1 in our 



definition of g, equation 4.46 



We take the gyrokinetic equation 4.7 for g and g' and multiply by (p' and 



ip respectively, then ensemble average. Next we integrate over velocity v. The 
result is 



(a - Tn 



(E) 



271 if' J vdv {{C[h])^)^ + 27ripJ vdv {{C'[h'])^,)^, + ^Tj + ^'T^ (4.54) 
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LfLf' and 



27r J vdv 5 {(f)^ 6\eSJ^ (with 



where we have defined $ 
^ {v)r = {'^')r' ~ (v^)r) which by homogeneity depends only on the increment 
i. Also by homogeneity, we have Fq^ = Fq^ = Fq^^$ where Fq^"* operates in i- 
space. Also, we have defined the i?-forcing JF^ = 27r J vdv {T)^. We can express 



the first term of equation 4.54 in terms of the rate of change of the invariant E 
and structure functions of y?. Using homogeneity we find the following relation 



[a 



- Fo)$ = - ^ (v - M(ro<^)) 



(4.55) 



where diV^Lp) = T'^ip' — Toip. We substitute this into equation 4.54 The structure 
function terms are zero in the stationary limit. With the additional assumption 



of small coUisionality, equation 4.54 becomes 



^TT c(^) , ^^r~^ , ^~Pr- 



(4.56) 



For £ = we have the global balance of E which describes the continual accumu- 
lation due to forcing: 



dE 
~dt 



(4.57) 



where Ee is the injection rate of E. Now recall that the forcing is assumed to 
be restricted to a scale ii as explained in the previous section on the forward 
cascade. For the inverse cascade, we are considering an inertial range such that 



i ^ ii- In this limit, the forcing term in equation 4.56 is zero and we may write 
the third order inertial range result for the inverse cascade of E: 



^An alternate equivalent expression for Sg^'is S'^^^ — Sip 2n J vdv {S^ESg) where 
Sip = — ip and the angle average with respect to r and r' is defined (/(R', R))^^. j.,-) = 
{27Tr'Jd^f{r'-p{^),r-p{^)). 
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(E) 



2e, 



(4.5J 



Using isotropy, we may also write 



(4.59) 



4.6.5 The CHM/Euler limit 



In the CHM or Euler limit (the result is applicable to both Euler and CHM since 



the electron response does not enter in the nonlinearity) , equation 4.58 reduces 
to 



(4.60) 



where we recall that vgo = z x Vy? is the E x B velocity without gyro-average. 



With some algebra, equation 4.60 may be manipulated to find agreement with 
the result from ( IGM02I1 for CHM turbulence and the result from (]Ber99p for two 
dimensional NS turbulence: 



(4.61) 



where u = — V£;o- 



The two gyrokinetic results of this section, equation 4.52 and equation 4.58 



hold for arbitrary spatial scales-i.e. they hold above, below and at the Larmor 



radius scale. (Note, however, that for equation 4.52 , we have assumed the velocity 
scale satisfies <^1.) In the the following sections, to derive approximate power- 
law spectra for the conserved quantities, we will consider the limit k ^ 1 and 
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p ^ 1. This limit corresponds to the regime where non-hnear phase mixing is 
dominant and gyrokinetics exhibits the phase-space cascade. 



4.7 The gyrokinetic turbulence cascade 
4.7.1 Self-similarity hypothesis 

For the forward cascade of free energy, scaling with respect to i is described by 



the property of S3, given by equation 4.53 As we have pointed out, self-similarity 



in is not predicted by the inertial range free energy-flux equation 4.53 Thus we 



will supplement equation 4.53 with an hypothesis as follows. Given fixed scales 



i and the difference 6g satisfies the dual self- similarity hypothesis 

Sg{Xi, XQ = X''<'6g{i, Q (4.62) 

by which we mean that Sg exhibits this scaling statistically-i.e. where it appears 
in an ensemble average. For instance it implies S2{Xi, Xiy)/ S2{i, iv) = X"^^^. The 
hypothesis may be justified from two perspectives. The first is the argument of 
decorrelation in velocity-space as has been described by (ISCDOSp and reiterated 



in section 4.5 A second point of view is based upon the fact that the collisionless 



gyrokinetic system is invariant to simultaneous scaling of R and v by the same 



factor (this is the symmetry 4.44d with 7 = 1). Thus, if we consider velocity 
increments ^t, -C 1 (that is, smaller than the thermal velocity, the characteristic 
scale of the Maxwellian), the statistics of the fiuctuations should depend weakly 



on the centred velocity v and obey a dual scaling in i and iy-i.e. equation |4.62 
As we will prove later, the scaling for 6g implies a self-similarity for the gyro- 
averaged E X B velocity difference ^v^;, for small i and i^: 
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5vij(A£, A4) = A^^5vij(£, Q (4.63) 



As a consequence of these scaling relations and equation 4.53 we see that the 



unknown constant in the expression for ^3 = dwE^g"^ must be linear in (if it is 
non-zero). Thus we write 

S,, = -2l{e(i + e'\Q) (4.64) 

where e' is an unknown constant. The absolute value on ^y is required, due to the 
fact that under the transformation (£, (i^) — > (— —(i^) we must have ^3 —S^. 
Finally, equation 4.64| implies that 



2hg + hE = l (4.65) 
4.7.2 The kinematics of gyrokinetic turbulence 

Before investigating spectral scaling, we must define the free energy spectral 
density. In this section we introduce a way of characterising velocity scales using 
a zeroth-order Hankel transform. The Hankel transform, being a two-dimensional 
Fourier transform with rotational symmetry, is also encountered in the theory of 
two-dimensional fluid turbulence. The difference is that in gyrokinetics, dynamics 
are exactly independent of the angle of the perpendicular velocity whereas in 
isotropic fluid turbulence, only statistically averaged quantities are independent 
of orientation in space. 

We now introduce some of the notation and identities that will be needed. 
Given a function g{v), its Hankel transform is 



137 



g{p) = / vdvJo{pv)g{v) (4.66) 







Inversion of the transformation is easily achieved using the orthogonahty of the 
zeroth-order Bessel functions, given by 



Jo{pix)Jo{p2x)xdx = S{pi- p2)/pi (4.67) 







Using this identity it is easy to see that the generahsed Parseval's theorem for 
two functions / and g is 

/ vdvf{v)g{v) = / pdpf{p)g{p) (4.68) 

JQ Jo 

An integral involving three Bessel functions is encountered in deriving the mode 
coupling due to the nonlinearity. This integral is 



/ vdvJo{piv)Jo{p2v)Jo{p3v) = K{pi,p2,P3) (4.69) 
Jo 

where if pi, p2 and p^ form the sides of a triangle, then K = l/27rA where A is 
the area of that triangle; if pi, p2 and ps do not form the sides of a triangle, then 
K = 0. Applying this identity, we can obtain the mode coupling equation for 
Hankel-Fourier modes. Neglecting collisions, the gyrokinetic equation, equation 



4.7 in Hankel-Fourier space is 



dg{k,p) 
dt 



^ j d'k'f3{k') (k - k') • (z X k') J qdq K{k',p, q) g{k', k')g{k - k' , q) (4.70) 
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4.7.3 The forward cascade of free energy Wg 

Returning to the problem at hand, the Hankel- Fourier transform of the distribu- 
tion function g is 

g{k,p) = ^ [ d^R [ vdvMpv)e-'''-^g{R,v) (4.71) 
Jr Jo 

Now we transform and inverse-transform g and use homogeneity to express the 
conserved quantity g"^ in terms of the correlation function G{£,v',v): 



g^ = J kdkpdpidiv'dv'Jo{k£)Jo{pv')Jo{pv)G{£,v',v) (4.72) 

where we have used the identities J jQ{ax)jQ{bx)xdx = 6{b—a)/b and J c/^ke*^'"^ = 
(27r)^5(r). Also, isotropy has been used to reduce the Fourier-transform in vector 
position increment £ to an Hankel transform in scalar increment i. Now recall 
the mean (ensemble-averaged) generalised free energy is defined 



Wg = 2TT I vdv— (4.73) 



From equation 4.72 it follows that 



Wg = J dkdpWg{k,p) (4.74) 



where 



Wg{k,p) = npk j idivdvv'dv'Jo{ki)Jo{pv')Jo{pv)G{i,v',v) (4.75) 
The spectral density can be expressed in terms of the transform function: 
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Wg{k,p) = npk g{k,pY (4.76) 



Notice that equation 4.75 is not invertible as with the case of the Weiner-Khinchin 
formula found in fluid turbulence-i.e. the correlation function G{i,v',v) cannot 
be recovered from Wg{k,p). Because of the lack of homogeneity in v, we have 
integrated over two dimensions v and v' to define scales in velocity-space via a 
single wavenumber p. 

4.7.4 Spectral scaling laws 

Now to establish the relationship between the scaling of the spectrum Wg{k,p) 
and that of the structure function 5*2 = Sg'^ at small scales, we will first change 
velocity variables to the increment variable iy = v' — v and the centred velocity 
V = {v + v')/2. 

Wg{k,p) = Tipk j idii^v^ - il)diydvJo{ki)Q{p,iy,v)G{i,i^,v) (4.77) 

where 



Q{p, 4, v) = Jo[p{2v + Q/2]Jo\p{2v - Q/2] (4.78) 

4.7.5 Scales smaller than the Larmor radius: A; ^ 1 

Until now, we have not made any assumption about the size of k or p. Now we 
now consider the phase- mixing range k ^ 1 and p ^ 1. First, we make use 
the fact that the correlation function G{iy) must be peaked around small to 
approximate that the integral is dominated by -C v. We can also take pv ^ 1 
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and use the large-argument approximation of the Bessel function. In this hmit, 
Q{p,£v,v) can be expressed 



Q{p, iv, v) ^ {sm{2pv) + cos(K.)) (4.79) 
Anpv 



Substituting this expression into equation 4.77 we see that the term proportional 



to sin(2pf ) oscillates rapidly and does not contribute in comparison, so we neglect 



this. Using G = {g'^ + {g'Y — 5'2)/2 we can express Wg{k,p) in terms of 5*2 for 
non-zero k and p: 

Wg{k,p) ^di vdv d£Mk() cos(K^,)^2 (4.80) 

Then self-similarity of 6g implies that Wg{k,p) must satisfy the scaling 

Wg{Xk,Xp) = X''Wg{k,p) (4.81) 

with 



z/ = -2 - 2hg (4.82) 
(Note that to arrive at this scaling, we have assumed that the bounds of the 



integration of in equation 4.80 are not important since we are concerned with 
the small iy behaviour of the integral.) To determine the scaling index u we will 
now find its relationship to He, the scaling index of the gyro-averaged E x B 



velocity structure function ^v^;. We will write S2E = I'^v^;^ in terms of the 
spectral free energy density Wg{k,p) and take the limit <^ v as above. We 
recall quasi-neutrality 
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0{k)=(3{k) J vdvJo{kv)K{v)g{k,v) (4.83) 

where 

m = (4.84) 
a - ro(fc^) 

and we recall that fo(/i;^) = e~*''^/o(/i;^). We momentarily assume the special case 
k{v) = 1. We must stress that this assumption will not, in the end, affect the 
generality of the results. After establishing the spectral scaling in Wg for the case 
K = 1, the scaling index of S2E is fixed, and one may refer to the general result of 



equation 4.64 to establish the general scaling of 5*2 and thus the general Wg{k,p) 
spectral scaling laws. The general-K scaling laws are found to be the same as 
those for to the case k = 1. Continuing, with k = 1, we may write 

m=m9i^,k) (4.85) 



which allows us to express S2E = l^v^;^ in the form 



S2E = J dkH{k,i,v,v')Wg{k,k) (4.86) 



where 



H{k,e,v,v') 



kf3\k) 



TT 



[Joikv') + Jl{kv) - 2Jo{kv')Jo{kv)Jo{ke)] (4.87) 



Expanding this in the limit kv ^ 1 and <^ v, as done for Q above, only part 



of H contributes significantly to the integral of equation 4.86 and we may take 
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H{k,i,i,,v) ^ [l-cos(H,)Jo(H)] 



(4. 



Substituting this into equation 4.86, we have 



2E 



dk [1 - cos{kQJo{ki)] Wg{k, k) 



(4.89) 



From equation 4.89 we may determine the scahng index in terms of z/, the 



scahng index of Wg-. 



2h, 



(4.90) 



Now we may combine our results for the high-/c hmit to obtain unique scahng 



indices for 6g, 6\e and Wg{k,p). From equations 4.65, 4.82 and 4.90 we have 



hg = 1/6 

= 2/3 
z/ = -7/3 



(4.91a) 
(4.91b) 
(4.91c) 



And, therefore, taking p = k we have the power law 



Wg{k,k) oc A;-^/3 



(4.92) 



which implies the large-/c limit of the second invariant E (which is simply the y?^ 
spectrum) scales like 



ak-, 



E{k) ^ aE^(k) ^ -^\0{k)\^ ~ k-^Wg{k, k) oc k-^^'^ (4.93) 
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This scaling for the electrostatic spectrum E^{k) is in agreement with the previous 



predictions of (1SCD07P (see section 4.5). 



4.7.6 The two-dimensional spectrum Wg{k,p) 

We have already seen that the two-dimensional spectrum Wg{k,p) obeys the 



scaling law given by equation 4.81, with u = —7/3. To obtain this scaling, we 
have assumed k ^ 1 and p ^ 1 but nothing about the relative sizes of k and 
p. In the subsidiary limits p ^ k and k p, we can show that the spectrum 
Wg{k,p) has a power law separately in k and p. The results we obtain are 

^-2^-1/3^ for k p 
Wg{k,p)^{ (4.94) 

p^^A;"^/^, for k <^p 

We will present the argument for p <^ k and omit the case k <^ p which is 
analogous. In the limit p <^ k, the scaling of the spectrum Wg{k,p) is determined 
by the scaling of the 5*2 in the limit that i iy. Thus we must determine the 
behaviour of the structure functions in this limit. Since we now have a power 



law for Wg{k,k), it is apparent from equations 4.89 and 4.92 that S2E = I'^v^;^ 
may be computed analytically as a function of i and iy. For general i and 
the structure function S2E is evaluated in terms of a hypergeometric function. 
However, here we are only interested in the limits i i^. Thus if we take 
Wg{k,k) = Wgo k~'^/'^ and expand the integral to first order in we obtain 
the result 

S2E ~ CdT + C2(i-^'^f (4.95) 
where Ci and C2 are constants which depend on WgQ and v. Next, we write 
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5*2 ~ \S3\/\/S2E and obtain to first order 



S, 



^1 



e t , 



(4.96) 



When we substitute this into the expression for Wg{k,p), equation 4.80, the 



leading order term produces a delta function after integration over i which will 
be zero for non-zero k. Since we are assuming k and p are both much greater 
than one, this term is obviously zero. Thus Wg{k,p) is 



Wgik,p) 



ke 



^1 



m vdv di^Joiki) cos{pQ (i-'^^^i) 



(4.97) 



It is easy to see that equation 4.97 implies the power law Wg{k,p) oc p ^^^k ^. 
Repeating this analysis for the limit k <^ p completes the demonstration of 



equation 4.94 



4.7.7 The one-dimensional spectra Wg{k) and Wg{p) 



By integrating equation 4.77 over p we obtain 



Wg{k) = I dpWg{k,p) = -7r| / idivdvJoiki)S2ii,v,v) (4.98) 



And integrating over k we obtain 



Wg{p) = j dkWg{k,p) = -\j vdvd£^cos{pQS2ii = 0,iv,v) (4.99) 
The scaling of 5*2 for zero gyro-position increment i can be obtained from equation 



4.95 The scaling of 5*2 for zero velocity increment iy is obtained analogously. We 



find the following power laws: 
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Wg{k) oc k-^^^ (4.100) 

and 

Wg{p) OC p-^/^ (4.101) 

The power law for Wg{k) has been predicted by (ISCD07P while the power law for 
Wg{p) is a new prediction. As a consistency check one may integrate Wg{k,p) di- 



rectly, using the asymptotic forms from equation 4.94, to obtain roughly Wg{p) 



Jq dkk ^/^p ^ + dkk "^p ~ p ^Z^. This works likewise for VI/g(/i;). 
4.7.8 The inverse cascade of E 

We now turn to the scaling of the spectra for the high-/c inverse cascade range, 



1 -C A; -C /^j. We recall the result of equation 4.59: S^^'^ = e^i (which is valid for 
= 0), where Sif ^ = 27r J vdv 6 {ip)-^6\ e^Q- 

We must take care in deducing the scaling index hg of the increment 5g from 



4.5 



we 



this expression for . From the phenomenological arguments of section 
expect that the velocity integration of 6g should, roughly speaking, introduce a 
factor of Thus we cannot, for instance, argue that the scaling indices satisfy 
hg + hE + hu) = 1 from the linearity of S-^ -as was done to obtain equation 



4.65 



i'g~r ii'E~r ii'{(f)) = J- uom Liie imeainy oi 

We will not use the random walk argument in this section. However, we 
will make the additional assumption of local transfer in fc-space (implicit in the 
phenomenology) to obtain an estimation of the effect of the velocity integration 
on scaling of S^^K In /c-space, 6g, 5\e and 5 {^p}^ are expressed 
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1 

1 

2^ 



g-iks-R' _ g-ika-R 



Jo{k2v)(p{k2) 

Joihv)izxk,)^{ks) (4.102) 



If we substitute these expressions this into the quantity 5 {(^)^SvESg, locahty 
in fc-space imphes that we may approximate Joik^v) ~ Jo{kiv) and Jo{k2v) ~ 
Jo{kiv) and transfer the two Bessel functions to the integral over ki. This leads 
to the approximation 



recalling that the E x B velocity (without gyro-average) is veo = z x V<y9. The 
scaling of 6ip and ^v^o may be written in terms of the index u by following the 
methods leading to the scaling index of S2E, h^. It is found that 



S^^Eom = A-^/2-i5vso(^) 



(4.104a) 
(4.104b) 



Now we determine the scaling of the remaining velocity integral in equation 4.103 



To do this we expand in the Hankel transform of 6g and perform the integration 



over velocity using equation 4.69 Then we square the quantity and ensemble 
average, applying homogeneity and isotropy. We obtain 
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Stt [ dkpdp K\k,k,p)[l- Jo{ki)]Wg{k,p) (4.105) 



where the function K is defined by equation 4.69 -it arises from the integral of 



three Bessel functions. Also, we have made the approximation that g{p)g{p') 



2TTg(p)g(p')5{p — p'). This approximation is analogous to the exact result in 



position-space due to homogeneity: g{k')g{\<:) = 6(k' + k)^(k')5f(k). The Hankel 
transform velocity-space analogue may be provec|^by taking the limit k ^ 1 and 
assuming that the correlation function G{iy) is highly peaked around iv = and 
varies smoothly in v (as we have been assuming in this chapter). 



4.59, 


4.103 


4.104 


and 


4.105 



iy = -l 



which implies, using equation 4.82, that 



(4.106) 



h, = -1/2 



(4.107) 



leading to the spectra 



E{k) oc k-^ (4.108a) 
Wg{k,k)ock~^ (4.108b) 
Wg{k) oc k° (4.108c) 



1 /-l/a 

'The identity d{p' — p) = lim — / cos[v{p' — p)]dv is useful here 
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4.7.9 Some comments on the two dimensional spectrum Wg{k,p) 



From the scaling relations obtained in this section we may write 6 {(p)-^6\'E6g oc 
However this result is specific to £y = 0. For ^ 0, the self- similarity hy- 



pothesis 4.62 can only give 5 {(p)-^ dwE^g oc /(£, such that the function /(£, (i^) 
obeys the scaling fi\£,XQ = \^/^f{i,Q. Without more specific knowledge of 
the function /, it is not possible to determine power law spectra for Wg{k,p) in 
the limits k ^ p and k <^ p. However, in the limit that i = and £^ 7^ we 
must have / oc il^"^, implying that 

Wg{p)ocp^ (4.109) 

4.8 Conclusion 
4.8.1 Summary 

In this chapter we have considered forced two dimensional gyrokinetics as a sim- 
ple paradigm for kinetic plasma turbulence. We have explored how the nonlin- 
ear phase-mixing mechanism gives rise to a phase-space cascade for scales much 
smaller than the Larmor radius. We have investigated the relationship between 
the inertial range cascades in fiuid theories (HM and two-dimensional Euler tur- 
bulence) and gyrokinetics, and found a simple relationship between the fiuid 
invariants and the kinetic invariants, given by equation |4.30[ The gyrokinetic 
free energy invariant Wg is cascaded to fine scales and is ultimately dissipated by 
collisions. Concurrently, there is a gyrokinetic invariant E which my in principle 
cascade inversely from very fine scales (e.g. the electron Larmor radius) to long 
wavelength regimes of the CHM/Euler limits, much larger that the ion Larmor 
radius. 
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In addition to a phenomenological derivation of the spectral scaling laws based 
on the previous work of ( ]SCD07p . we have given derivations of exact third-order 
laws for the forward cascade range and for the inverse cascade range. These 
derivations are based upon general considerations of symmetries of the gyrokinetic 
equation and an additional assumption of a two scale dependence in velocity- 
space. 

Combining the exact third order results with a dual self-similarity hypothesis 
and the two-scale assumption in velocity-space, we are able to reproduce the 
phenomenological scaling laws for the both the forward cascade range and the 
inverse cascade range. Since the methods used are somewhat different, this can 
be seen as an independent confirmation of these scaling laws. 

For the forward cascade range, the dual self-similarity hypothesis is also used 
to obtain novel predictions for the velocity-space spectrum. To describe the 
velocity spectrum, we introduce the use of a zeroth order Hankel transformation. 
This novel spectral treatment of perpendicular velocity-space may in general, be 
found useful for theoretical and numerical applications in gyrokinetics. 

4.8.2 A note on forcing and universality 

In general, turbulence may be generated by linear instabilities which are induced 
by a mean-scale free energy gradient. In this case, additional (linear) terms in 
the gyrokinetic equation must be included. Although the present work has left 
the exact form of the forcing unspecified, we note that a simple phenomenological 
calculation shows that the additional terms do not invalidate the inertial range 
assumption at fine scales - the details of these arguments are left for future work. 
Thus it is possible that these scalings may be a universal feature of magnetised 
plasma turbulence. 
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Encouragingly, the scaling for E^f,{k) in the forward cascade range has already 
been independently observed in a fusion- relevant simulation by (IGJOSap . Also, 
the simulations found in fITBCOSP demonstrate strong confirmation of the scalings 
predicted by (1SCD07I1 and also demonstrates consistency with the velocity-space 



prediction for Wg{p), given by equation 4.101, Numerical investigation of other 
novel predictions of this chapter, including the inverse cascade scaling laws and 



the detailed phase-space spectral predictions of section 4.7.6 are planned. 



4.8.3 Future work 

There are some issues not addressed by this work which should be given high pri- 
ority. First, a more complete investigation of the gyrokinetic turbulence cascade 
would include physical content such as electromagnetic affects and the correct lin- 
ear instability drives. Electromagnetism introduces a fluctuating magnetic field 
into the nonlinearity which couples to the distribution function through Ampere's 
law and will enrich the problem substantially. 

The present work requires significant extension to be applied to the realistic 
tokamak scenarios. As a general point, it should be noted that the macroscopic 
properties of a fusion plasma are most strongly influenced by scales at the forcing 
range or larger. This has meant that inertial range ideas are generally not a focus 
of modern magnetically confined fusion research. However, some more immediate 
applications of this work to fusion include simulation benchmark and diagnostics, 
or the general use of the v± Hankel-transform spectral treatment in gyrokinetic 
linear and nonlinear problems. 

Inertial range understanding gyrokinetic turbulence may also enter in the 
problem of transport. As turbulence modelling in fusion ultimately aims to pro- 
duce accurate quantitative prediction of transport, the inertial range contribution 
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to transport will likely need to be included as the field evolves. 

Three dimensional dynamics must also be investigated. First, the assumption 



that the parallel streaming term (see equation 4.2) may be neglected requires 
careful justification. One instance where this term does play a role is in the 
Alfven wave cascade-which can be treated within gyrokinetic theory. In this case, 
it is believed that a critical balance condition is satisfied where the generation of 
finer scales in fcy ensures that the parallel term remains important at all scales in 
the cascade forward. On the other hand, a conventional assumption in tokamak 
turbulence is that one may take the size of to be correspond to the distance 
along a field line from the stable side to the unstable side of the toroidal magnetic 
surface. In that case, the nonlinearity will dominate over this term for sufficiently 
small (perpendicular) scales. 

It is known that the parallel streaming term also gives rise to a linear phase- 
mixing, which generates fine scales in parallel velocity-space (fy) dependence of 
the distribution function. Nonlinear simulations and theoretical work by (IWS04p 
demonstrate how this parallel phase-mixing plays a role in the statistical steady 
state of the entropy cascade. Their work describes a steady state where parallel 
phase mixing is treated as a passive scalar phenomenon and exhibits an inertial 
range with a power law spectrum. However, this work is done employing a v±- 
integrated version of the gyrokinetic equation, which suppresses perpendicular 
phase-mixing by assuming a fixed Maxwellian dependence in f ^-space. Thus our 
work is complimentary to the theoretical work by ( ]WS04p . although both works 
neglect inclusion of the third dimension in position-space. 

In general, the present capabilities of numerical simulations allow a detailed 
investigation of (perpendicular) nonlinear phase mixing or (parallel) linear phase- 
mixing, but not both. Thus a theoretical investigation of the full five-dimensional 
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theory (that is, three position dimensions and two velocity dimensions) currently 
has the opportunity to pave the way to a more complete picture of the turbulent 
steady state in the gyrokinetic turbulence cascade. 
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CHAPTER 5 



Conclusion 

5.1 Highlights 

In this thesis we have seen the apphcation of varying approaches to modehng 
and understanding the fully developed state of a turbulent magnetized plasma. 
There has been a common motivation, in magnetically confined fusion; and there 
has been a common theme, in the multiple scales at play. In chapters [2| |3] 
and |4] we have seen this theme in the multiple scales approach to gyrokinetic 
transport theory, in disparate time-scales in the saturation of linear instabilities 
by secondaries, and in the treatment of velocity scales in the phase-space cascade 
with a two-scale assumption. 

Chapter |2] (and the appendix |A]) has introduced a synthesis of gyrokinetic and 
(neo) classical transport theories based on the method of multiple scales. Chapter 
|3] has explored secondary instability theory as a simple model of an important 
nonlinear process in fully developed turbulence. We have argued how this the- 
ory can be applied to mixing length phenomenology to shed light on important 
broader issues of tokamak physics. Finally in chapter |4] we have developed a 
framework for understanding the turbulent generation of phase-space structures 
and derived the inertial range scaling laws, following the methodology from neu- 
tral fluid turbulence. 
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5.2 Final note 



Although each of the works in this thesis have produced novel results, they are 
intermediate steps in solving larger and more difficult problems - and, naturally, 
there are important questions and opportunities which remain. The chapters 
of this thesis have already detailed some of the future work that is needed. In 
the context of the combined work of this thesis, there are some general (and 
somewhat speculative) remarks that can be made. 

Statistical theories of fully developed strong turbulence in the energy-containing 
range have been unable to capture the features that were subsequently observed 
in nonlinear simulations such as the levels of anisotropy and intensity and, ulti- 
mately, the thermal fluxes. These analytical theories are based upon plausible 
assumptions about the nature of the fully developed state. 

For instance, a typical assumption is that, at each scale, the linear growth rate 
of unstable modes should be balanced against the nonlinear decorrelation time 
determined generically by the strength of the nonlinear term. (Or more generally, 
one may assume that the nonlinear decorrelation time may be determined as 
a function of the linear growth rate spectrum - see (ICGC87p .) However, for 
important cases in tokamak turbulence, it has turned out that this plausible 
assumption misses the mark. For instance, although the linear theories of ETG 
and ITG turbulence are isomorphic, the saturated turbulent state in the energy- 
containing range are qualitatively very different, as discussed in chapter |3| 

Secondary instability theory has shed light on this issue. For the case of 
ETG turbulence it can be shown that, for sufficiently small wavenumbers, the 
magnitude of the nonlinearity required to destabilize a fast secondary becomes 
much larger than the level needed for a basic balance with the linear drive terms. 
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These weak secondaries have been used to successfully predict the shape of the 
saturated spectrum (see ( ]JD02p ) for an important part of the energy- containing 
range. That is, it turns out that the shape of the ETG saturated spectrum is 
largely determined by the balance between secondary and primary growth rates. 

As numerical simulations and focused analytical models build our intuition 
for the basic physical processes which govern the different turbulent ranges, there 
will be an opportunity to revisit the problem of a more complete and fundamental 
theoretical description. Such a description must contain essential ingredients of 
the linear and nonlinear physics which leads to such things as self organization (or 
spectral condensation) and observed saturation amplitudes. This theory would 
be a powerful tool in creating a deeper understanding of fusion confinement. 
It could lead to a phenomenology of magnetized plasma turbulence with broad 
applicability and, in general, a mathematically and computationally simplified 
framework to enable efficient exploration of the vast parameter space for the 
operation of fusion devices. Finally it may be hoped that this would be a path to 
the discovery of ways to tame turbulent plasma transport and realize the dream 
of small-sized and inexpensive fusion reactors. 
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APPENDIX A 



Appendix: Gyrokinetics transport in toroidal 

geometry 

A.l Introduction 

The appendix is devoted to a lengthy derivation of transport and entropy balance 
for axi-symmetric toroidal geometry. The procedure followed here is quite similar 
to that in chapter |2| so there is much repetition. The appendix has not been 
produced in a form designed for publication in a journal (the work is still in 
progress), so it should be viewed as a collection of notes, though having several 
concrete results. 

Others who have combined neoclassical and anomalous transport are (lSha88| 
IBalQOt ISOH96al) . The work which is probably closest to ours is flSOH96aP . They 
also derive the transport equations within the gyrokinetic ordering hierarchy, but 
employ an ensemble average to separate the mean-behavior, instead of the patch 
and time averages we use based on scale separations. Also, the final form of heat 
transport we have provided is significantly manipulated as to compare closely 
with entropy balance, and so that the heating term is in a positive-definite form 
suitable for application to numerical simulations - see (IBarOSp . 
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A. 2 Equilibrium field geometry 
A. 2.1 Axi-symmetry 

We assume toroidal geometry and make use of axi-symmetry when needed. How- 
ever, we avoid explicitly using magnetic coordinates when possible. We make 
the common choice of labeling magnetic surfaces by the flux variable which is 
related to the poloidal flux by = ■0p/27r. The equilibrium magnetic field can 
then be expressed 



where (p is the toroidal angle and is a flux surface quantity which is propor- 
tional to the total poloidal current. 

A. 2. 2 Magnetic coordinates 

There are many suitable choices for flux coordinates. Our choice is the radial 
coordinate ip, the poloidal angle coordinate 6 and the toroidal angle (j). There is 
also some freedom in choosing the definition of 6. Here we define 



where Ip is the length along the flux surface as one integrates in the poloidal 
direction in a constant plane. This choice of 6 has the feature that field lines 
appear "straight" in the (f)-9 plane. In particular, a field line satisfies the equation 
^ — q, where ^ is a constant on the flux surface called the safety factor. (The 
safety factor is the ratio of the number of toroidal transits to the number of 
poloidal transits that a single field line makes as it traces out the fiux surface.) 



(A.l) 




(A.2) 
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These coordinates provide another simple way to express Bq 



Bo = VV^ X (0 - qe) (A.3) 
As a final note, the equilibrium field direction will be denoted by 



bo = 1^ (A.4) 



A. 2. 3 Annulus volume average 

The annulus volume average is defined over an interval Atp containing the fiux 
surface with label ip- In particular it is defined 

where 

/ d^rA= / / / — ^-^——A(ij',9,^) A.6 

Jav Jo Jo VV^'xV^-V0 '"^^ ^ ' 

The annulus volume average is inspired and closely related to the fiux surface 
average from neoclassical transport theory. In fact, when applied to equilibrium 
quantities, it is equivalent to the fiux surface average since the shell can be consid- 
ered to have infinitesimal thickness. Indeed, the explicit form of the fiux surface 
average in Ref. (lHH76p . Eqn. 2.55, is identical to the annulus volume definition 
in the limit of infinitesimal volume AV. We will make use of this fact in adapting 
results from neoclassical theory, Refs. (1HH76[ IBer74p . For the turbulent terms, 
the annulus has finite thickness and essentially combines a fiux surface average 



with the patch volume average defined by Eqn. 2.15 The annulus volume average 
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of the divergence of a vector function F is evaluated in a slightly different manner 
to the flux surface average 



1 



(V ■ F)^^ = 



Ay dtp J iwi ^ ^ 
^^(Ay(F-V^)^^) = 

dI^^ (F ■ V^)^^ (A.7) 

where we have introduced the operator d[^^ defined D^^^A = AV^^d/dip{AV A). 
I will also make use of the analogous result for the flux surface average (Ref. (lHH76p . 
Eqn. 2.58) which I copy here: 



where V is the total toroidal volume contained within the flux surface labeled ip, 
and the we define D^A = dip/dV d/dip{dV/dipA). 



A. 3 Proving Fq is a flux surface quantity 



In Sec. |2.6.2| we used the assumption that the equilibrium distribution function 
does not vary in the direction of the equilibrium magnetic field, bo-V-Fo = 0. This 
can be proven for a system with closed flux surface geometry (e.g. axi-symmetric 



toroidal geometry). First we recall the 0{1) Fokker-Planck equation, Eqn. 2.31 
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C(Fo,Fo) (A.9) 



We use Boltzman's H-theorem to show that Fq must be MaxwelUan. We multiply 



Eqn (2.31) by 1 + logFo, integrate over velocity space, but this time integrate 
over annulus volume defined in Sec. IA.2.31 The troublesome term is evaluated as 
follows: 



(A.IO) 



where we have used the divergence theorem to evaluate the volume integral as 
a surface integral and used the fact that the surface vector is perpendicular to 
the equilibrium magnetic field so that d'^s ■ v = rf^s ■ v^. This is obviously odd 
in "i? since Fq is independent of So, the only term remaining from the 0{1) 
equation is the coUisional generation of entropy per unit volume in the annulus 
volume. So we have 

^y"(i3vC(Fo,Fo)lnFo^ =0 (A.ll) 

As before we have that Fq is a Maxwellian. Finally, we return to the 0{1) 
equation and apply the patch average and average over i) leaving the desired 
result: 

bo-VFo = (A. 12) 
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Since we are considering a closed flux surface geometry (and ergodic field lines), 
the lack of variation along the magnetic field implies that Fq is constant on a 
fixed flux surface. I.e. in space, Fq is a function only of the radial coordinate ■0: 



(A.13) 



A. 4 0{e^): Transport equations 

Here we calculate the particle transport and heat transport equations in Gy- 
rokinetics. To obtain the long time dependence of Uq and Tq we consider the 
moment equations of the full un-gyro-averaged kinetic equation for any species 
(5 = Csrifs, fr)+Cssifs, fs) where Csrifs, fr) IS coUisious of species s on species 
r and Css{fs, fs) is like particle collisions). This avoids having to write out all 
the 0{e^) terms from the expanded equations. 

A.4.1 Pcirticle transport 

To obtain the time evolution of the density, we integrate the Fokker-Planck equa- 
tion over velocity and average over the annulus volume. 



The integral of the collisions over velocity will be zero to conserve particles. The 
(E + V X B) • 1^ terms will be zero because they are a perfect divergence in 
velocity space. What is left is the continuity equation. 




dt 



+ v-V/. + ^(E + vxB)--^ 




(A. 14) 
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dt 



+ V ■ (v/. 



Ay 



The flux term is evaluated using Eqn. |A.7 



(A.15) 



(A.16) 



We can transform this expression by considering the following 



V X Bo 



df 

dv 



AV 



V X BofR 



AV 



d\ (v ■ Vipf) 



AV 



(A.17) 



where we first integrated by parts, and then used Bq = Vip x V0 + VipF^ip)- At 
this point, we can again substitute the Fokker-Plank equation in for v x Bq ■ 
Thus, we can write 



/ rf^v (V ■ (v/,))^^ = / d\ ((v ■ V^/))^^ 



<V^^v(fv.0[^vxBo-f] 



AV 



/ ^'v (^v ■ [f + V ■ V/. + ^(E + V X <5B) • f - CU; /)] /A. 18) 



We can rewrite Eqn (A.15) as 



'AV 



Id'- it). 

^(A) J rf3^ (^11^ + V ■ V/. + ^(E + V X 5B) ■ - C{f, /)] )^^A.19) 



At this point, we would use our solution for the distribution function given by 



Eqn (2.44) and substitute it into Eqn (A. 19). Most terms are zero or negligible 
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by averaging and ordering arguments. We leave the work for Sec. A.7 The resuh 
is 



/ dnps \ 
\ dt / 



AV 



Df^ J d^^r (r4> ■ + Vx^k + fC{p- VF^ - h^c, F^) ) (A.20) 



Ay 



where Fm is the MaxweUian. The first term is the Ware pinch. The second term 
is accounts for the turbulent or anomalous transport driven by the fluctuating 
fields and the associated drift velocities. The coUisional term contains classical 
and neoclassical particle transport. 

A. 4. 2 Electron particle transport 

The above equation and analysis in Sec. A.7 is independent of species. The 



density transport for electrons is simplified by integrating the electron gyrokinetic 
and neoclassical equations over velocity. The result obtained is that hgk and hnc 
are both odd v\\. The resultant electron density transport equation reads 



' dn 



Oe 



dt 



D 



(A) 



Ay 



(A.21) 



Ay 



Note if we assume the Boltzmann response for electrons, the anomalous term 
(the second term) will integrate away which implies there will be no anomalous 
particle transport of electrons for ITG-driven turbulence insofar as the Boltzmann 
response holds. 
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A. 4. 3 Heat transport 

To calculate heat transport, we multiply the Fokker-Planck equation by |mv^, 
integrate over velocity, and average. Following the same methodology as with 
the particle transport. 



mv 



|+,.V/+l(E + vxB).|^ 

ot m av 



AV 



(A.22) 



AV 



The second term in Eqn (A.22) can be rewritten in the same manner as Eqn (A. 17). 



/rf»v=f!(v.V/>,, = C{-'/d»v!f!(v.V^*/> 



Ay 



(A) 



D 



mv'^R 



AV 

df 



— IV X Bo; ■ 
m av 



(A.23) 



Ay 



The last line was obtained using (v x bg) ■ ^ = 0. Reworking of Eqn (A.22) 
yields 



mv"^ df \ 



2 dt / 



AV 



2^ 



+ / d^v 



mv^R 



^(v X Bo) . ^ 

m av 



gE . V/ + /) 



/ Ay 
(A.24) 



Ay 



Like the particle transport analysis, using Eqn (2.44) in Eqn (A.24) is left for 



Sec. (A. 8). The result we obtain is 
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'3d{nosTos) 
2 



dt 



-D 



(A) 



FmV d) ■ — — 

2 ^ dt 



AV 



D 



(A) 



^— (Vx • (t>)hgk + ^^:^v • 0C(/i„c - P ■ VF^, F„) 



/ 



dt 



h 



gk 



2q 

dAo dAo, 



Ay 



dt 



dt 



Ay 



(A.25) 



The first term will give us the temperature evolution. The second through fourth 
terms are heat convected by the particle transport. The ^hgk term is gyrokinetic 
heating. The second and third terms on hne three are heating from the Ware 
pinch. The last term is heat exchange between species. 



A.4.4 Electron temperature transport 

Because hg^f. and hnce both are odd in ^n, the above heat transport equation is 
again simplified for electrons. 



3 dinoeToe) 
2 



dt 



D. 



(A) 



-D. 

AV 



(A) 



rripR 



FmeV (p ■ 

dt 



AV 



AV 



dip dAo ■ 

-Q^hgke-V -Q^{p-VF^e) 



noeiyl'{Ti-T,) (A.26) 



Ay 



A. 5 Closure of the moment equations 



In this section we describe the necessary steps to obtain closure and complete 
the transport time step to obtain the equilibrium evolution. At this point we 
have the equations for turbulence (the electron and ion gyro-kinetic equations), 
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and the equations for transport (heat and density). The prescription for solving 
this system is to evolve the gyro-kinetic distribution function over many turbulent 
time steps and use the averaged solution as input for the transport equations. 
In other words the turbulence should evolve to generate an average steady-state 
turbulent transport (if, however, the transport grew unbounded and did not 
reach this steady-state, the gyro-kinetic approximation would become invalid al- 
together). In principle one would then proceed to evolve the density and pressure 
(temperature) profiles a transport timestep to obtain the new equilibrium profiles. 
However, this step is not quite complete because we have left to solve for the time 
evolution of the equilibrium field geometry. We now show how the evolution of 
the equilibrium geometry is obtained to complete closure. 

A. 5.1 Flux surface motion 

The partial time derivatives such as dno/dt have been thus far implicitly taken at 
fixed r. However, since the temperature and density are flux surface quantities 
and the flux surfaces themselves move on the turbulent time scale, it can be 
useful to know how density and temperature of each species evolves on a fixed 
flux surface. The following argument can be found in Ref. (IBer74p . For a fixed 
point in space, the motion of flux surfaces must satisfy a continuity equation of 
the following form: 

^+v^-VV^ = (A.27) 

where is the flux surface velocity. Since motion parallel to surfaces is not 
physically meaningful we can write without loss of generality that 



168 



" 9^ |V^|2 



Thus we arrive at the following identity 



(A.28) 



Combining this result with Eqn. |A.8| we have 



(A.29) 



r/ ip 



(A.30) 



An expression for ^ follows directly from taking the inner produce of Vip with 
Faraday's Law for the equilibrium field —dBo/dt = V x Eoi (see also (1HH76P ) 



dip 
~dt 



where we define Eoi 



-R(j) ■ Eoi (A.31) 
-dAo/dt, the equilibrium inductive electric field. We now 



use Eqn. |A.30| and Eqn. |A.31| to rewrite the transport equations. Using the fact 



that density and temperature are flux surface quantities we write 



dt 



AV 



dnps 
dt 



Os 



dt 



- Dj, {nosR(j) ■ Eoi 



(A.32) 



and 



3 / dnpsTps 
2 \ dt 



AV 



3 ^ dnosTos 
2 



dt 



3 dnpsTps 
2 dt 



( 7T^Os^Os-R0 ■ Eoi 



(A.33) 



Substituting these expressions into the density (Eqn. A. 20) and heat (Eqn. A. 52) 



transport equations gives cancellation with the ware pinch terms: 
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Os 



dn, 
~dt 



/ d\ I R(t) 



mv 



^Xhgk + — C{p ■ VFm - hnc, Fm) 



(A.34) 



and 



3 dnpsT ^ 
2 dt 



Os 



d■'^r 



mv'^R 



hgkT 



m 



dF 

dijj 



hgkT 



q dip 



- T 



dt dip , ^ 



(A.35) 



Notice that we have adopted the flux surface average in place of the annulus 
volume average for all terms in the transport equations (including the turbulent 
terms). In fact, the formal smoothing procedure we have been employing thus far 
is more than needed. The toroidal symmetry of our system of equations ensures 
that the flux surface average will be sufficient as an effective ensemble average. In 
other words subsequent smoothing over intermediate time and space (A?/)) inter- 
vals are redundant and we may compactly express the final transport equations 
only in terms of the fiux surface average from neoclassical theory. Consider now 



the form of the equilibrium magnetic field given bye Eqn. A.l in Sec. A. 2.1 



Bo = VV- X V0 + /(^) V0 
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The equilibrium field is determined when one knows the spatial dependence of the 
poloidal flux function ip{R, Z) and also the toroidal field function It is a nat- 

ural and common choice to therefore solve for these quantities. This is essentially 
done using Faraday's Law and the Grad Shafranov Equation. The equations ob- 
tained, however, depend on the flux volume averaged quantity (Bo ■ Eoi)^^ which 



must be obtained by solving the neoclassical equations Eqn. 2.54 and Eqn. 2.56 



A. 5. 2 Evolution of the toroidal field function /(^) 

Our stated aim is to now solve for the time evolution of functions and 
ipiR, Z). We follow Bernstein (lBer74p and take the toroidal component of Fara- 



day's Law to obtain 

a/ 

Applying the flux surface average to this equation then yields 



^^i?-2 = V-(V0xEoi) (A.36) 



= /(Eoi ■ Bo - ^i?-^/) ) (A.37) 



Finally, the last term will drop out upon using Eqn. A. 30 



Which can be written in an alternate form by using the chain rule on the left 



hand side and another application of Eqn. A. 30 
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This equation gives a way to solve for /(V^) at the next time step. Combined with 
the Grad Shafranov equation and the heat transport equation (Eqn. A. 35) one 
can calculate ip{R, Z) at the next time step. The Grad Shafranov equation is 



R'V ■ (R-'V^P) = -I^ - i^oR'^ (A.40) 



where 



s 

is the total equilibrium pressure. However, we are not quite done because we 
have not yet shown how one obtains the quantity (Bo ■ Eoi)^. 

A. 5. 3 The neoclassical equation and obtaining (Bq • Eoi)^ 

The quantity (Bq ■ Eoi)^ is not known a priori. It depends on the inductive equi- 
librium field which is dependent on the evolution of the equilibrium magnetic 
field for which we are seeking a solution. The quantity (Bq ■ Eoi)^ can be deter- 
mined, however, by solving the neoclassical equation. One uses the solution to 
calculate the flux volume averaged parallel current which is known by Ampere's 
law in terms of spatial derivatives of the equilibrium magnetic current. We shall 
see that this gives an expression for (Bq ■ Eoi)^ in terms of the solution to the 
neoclassical equation and known quantities. We recall the neoclasical equation 
for ions (electrons are treated analogously) 



^;||bo ■ Vhnc - C{hnc) = -VD ■ VF^ + g^^^ijbo ■ Eoi (A.42) 



The approach (see (IHH76P and (IBer74p ) is to transform hnc so as to obtain an 
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equation involving (Bq ■ Eoi)^ and not the full Bq ■ Eoi. The transformation we 
will use is 



h + 



qFm 

To 



de' 



bo ■ we' 



bo ■ Eoi — Bq 



'Bo ■ Eoi)^ 



(A.43) 



where e is the poloidal magnetic coordinate defined in Sec. |A.2.2[ Note that the 



second term in the integral is determined by the requirement that hnc is peri- 
odic in e. Note also that the integral term is proportional to the maxwellian 
of the species. If one examens the transport equations it is apparent that this 
maxwellian part will yield zero since it in coUisional terms involving other maxwellians. 
Therefore, we will be able to use hnc in place of hnc wherever it appears. Due 
to axi-symmetry, we can write the first term of the neoclassical equation as 
?;||bo ■ V6'^^. The resulting equation for hnc is 



v\\ho ■ Vhnc + C{h 

r 



-vz,.VF„H-«„«^i^B„<?^^ (A.44) 



Due to linearity, we may further split up the solution into two pieces for the two 
driving terms. 



hnc = 9l (Bq ■ Eoi)^ + 92 

where gi and g2 separately satisfy the equations 



(A.45) 



qFm Bq 



vi\ho ■ Vgi - C{gi) = 
v\\ho ■ Vg2 - C{g2) = -vd • VF, 



To {BD^ 



(A.46) 



These equations must be solved to obtain hnc- The flux volume averaged parallel 
current can be expressed in terms of these solutions. 
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ehr 



(A.47) 



this is then equated to an expression for {J\\)^ obtained from the Grad-Shafranov 
Equation and the parallel component of Ampere's Law 



The result is an expression for (Bq ■ Eoi)^ 



dl_ 



(A.48) 



(Bo • Eoi) 



A. 5. 4 Closure 



dl 



TJ _i_ J_dP\ 

aV^O Bo 



(/c/3v v\\{qgii - egie)) ^ 



(A.49) 



The equations needed for closure of transport are Eqn. A.34 , Eqn. A.35 Eqn. A.38 
Eqn. |A.40 and Eqn. A.49 The density transport equation can be integrated a 
time step independently. The others, however are mutually dependent and must 
be solved simultaneously and consistently to obtain Piip), and ip{R, Z) 



at the next time step. The heat transport equation (Eqn. A.35) depends on 
{dip/dt)^ and the evolution equation for /(■?/') depends on {R^'^dtp /dt)^. The 
value of dip/dt obtained by evolving ip{R, Z,t) using the Grad Shafranov equa- 



tion must therefore be checked for consistency with the value used in Eqn. A.38 
and Eqn. |A.49 
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A. 6 Entropy balance 



The purpose of this section is to calculate the entropy production of the plasma. 
Standard theory predicts that the time change of the local entropy will be due 
to either convection or collisions. 



^)^^ = ^r^/ rf^v(V^-(v/ln/))^^- I d\{C{f,f)lnf)^y (A.50) 



By manipulating the left hand side of Eqn (A.50), we can insert the density and 



temperature evolution equations we have derived in Sec. A. 4 



/dS\ 



AV 

dn, 



3 / diyFrnJ'Tf'Ffy 



dt 



AV 



dt 



AV 



3 , ,m . 
lnno + -ln(- 



lnT + 1 



3 /^KTo) 



2T \ dt 



(A.5i; 



Our goal will be to show that the terms we have calculated for the density 
and temperature transport indeed produce only convection or entropy produc- 
tion through collisions. In order to achieve this result, we must manipulate 



Eqn (A. 25), the temperature transport equation. The details of the manipula- 



tion can be found in Sec. (A. 8.1). 
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3 d{nosTi 



Os) 



2 dt Uy 



dt 



Ay 



+ <)/.3v(!^(Vx-0)V.^ 



2 



; Pm) 

^1 I C\V 



hgkT 



d\ { -^R(t) ■ Vx 



{C{hgk)) 



R 



Ay 



q dip 
dAo dF„ 



Ay 



dt d^/j I 



(A.52) 



We can now substitute in from Eqns (A. 20) and (A.52) into Eqn (A. 51) 



, / Ay 
d\hnno + -H- 



lnT + 1 



dAo 

dt 



AV 



+ 



T 



+ / a V ( — R<p ■ \-Vxhgk - -gf^rn H ^ — C{Kc 



+ / d\ 



dip 
h 



{C{hnc, Fm)) 



{C{hgk,Fm)) 



AV 

nosys—^ 



Ay 



AyJ 



(A.53) 



At the moment, the last hne of Eqn (A.53) contains terms that are neither con 



vection or colhsional. What is left to do is pull out the -Oi^'' to act on the entire 
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second line, put the ^ inside the D^^^ of the second hne, and evaluate dFm/dip 
on the third line. We find 



Ay 



:^v/^(C(5/i,F„)) 



Ay 



(A.54) 



where we take 6fi = —p ■ VFm + hnc + hgk, noting that here the Boltzmann part 
does not contribute to collisional entropy production. This form is exactly what 



we expect from Eqn. |A.50[ The last two terms can be easily identified as the 
perturbative expansion of fC{f,f) with the final term giving the inter-species 
collisional entropy exchange. 



A. 7 Particle transport calculation details 



This section will demonstrate how to obtain Eqn (A.20). We begin with Eqn (A. 19) 



j3 / 9fs 



Ay 
mi?v ■ 



^ + v-V/. + ^(E + vx5B).§^-C(/,/) 
ot m av 



Ay 
(A.55) 



Our goal will be to use the solution for the distribution function given by Eqn (2.44 ) 
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and evaluate the terms in Eqn (A. 55). Plugging the solution for / into the left 
hand side of Eqn (A. 55), we obtain up to O(e^) 



'dfs 
. dt 



AV 



dnos\ ^ /■ ,3 /dh\ ^ /• ,3 /d5f2 
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We now have in our equation ^7^, which is the time evolution of the density. 



The remaining terms in Eqn |A.56 will drop out upon time averaging. The first 
term on the right hand side of Eqn. A.55| is treated analogously to Eqn. A.56 but 
the order is one smaller because of the preceding factors. Thus, it contributes 



nothing. The next term to calculate in Eqn. A. 19 is also zero as we will see. 



rf3^(!!!:?v-0(vV)(Fo + /i + 5/2) 

(1 / AV 



(A.57) 



The 5/2 can be rewritten as a perfect divergence to the order we are calculating 
and will drop in ordering after averaging it over space. The Fq term is odd in 
velocity space and will integrate away. The gyrokinetic part of the h term will 
space averages to zero, and all that is left is a neoclassical term. We rewrite this 
by defining a symmetric pressure tensor P. 
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(A.58) 



so that we have 



^ j d\ i^^R ■ (V ■ P))^^ = ^ j d\ ■ (P • - P : L^^^ (A.59) 

where L is the antisymmetric tensor V(0R) = 0R — R0. (Note: R is the radial 
unit vector in cylindrical coordinates.) The second term is identically zero being 
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the double dot product of a symmetric and an antisymmetric tensor. The first 
term is a perfect divergence and we again use the divergence theorem to evaluate 
it. Referring to the above expression for P we see that P ■ is everywhere parallel 
to the flux surface (perpendicular to the surface element vector) so the divergence 
theorem integral yields zero. The field terms in Eqn. |A.55 can be written as 



d^^ / —^r ■ 0^(E + V X 5B) ■ ^(Fo + h)\ (A.60) 
\ q m ^^r / 

When considering Fq terms, the —p-VFm term is shown to be too small in order 
since the preceding perturbed fields have a small patch average. Since oc vFm 
we worry only about the electric field for the Fm and the 6fi Boltzman correction 
terms and write 



^ q m o\ I 
I (9A 



AV 



dt 



AV 



Were we have used that -R0 ■ Vy? integrates to zero over the annulus. What is left 
will time average to be negligible. We now write E + v x 5B as — Vx — v ■ V5A 
and calculate the h term by integrating by parts in velocity, so that 

d 

AV 



{km) 

where the v • VA term is removed in two parts. First the f ||bo • V5A is written as 
a total divergence with negligible error and vanishes due to divergence theorem. 



Second, the ■ V5A is removed by employing Eqn. 2.41 after integration by 
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parts. (Note that what we are interested in is the gyrokinetic part of h-the term 
with hnc does not enter at this order, because the patch average of Vx and v ■ V A 



are smalL) The final term in Eqn. A. 55 simply involves plugging in values for the 
distribution function. 

' mR 



I 



1 / Al/ 



(A.63) 



Recall that we denote as the Maxwellian so as to distinguish it from the 



modified Maxwellian Fq defined in Eqn. 2.40 Collecting terms, and recalling 



the -DI^^ in front of the integral inside Eqns (lA.Gll), (|A.62l), and (|A.63l) from 



Eqn (A. 55) , we have the transport equation for particles. 
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A. 8 Heat transport calculation details 



In this section, our goal will is to substitute our solution for the distribution 



function into Eqn (A. 24). 
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We begin with the left hand side of Eqn (A. 65) will yield 
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(A.65) 
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I d\ (^Irn^r'^^ih + 5/2)^^^ (A66) 



180 



The other terms in the time derivative will be negligible after time average. The 
first term on the right hand side of Eqn (A. 65) will be treated very similar to the 
particle transport starting with Eqn (A. 18). 
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The significant difference between this equation and the particle transport term is 
obviously the inside the integral. Fortunately, all the terms that were negligible 
in the density transport will still be negligible as we will see. The time derivative 
will time average away. We can remove all parts of the next term except the 
neoclassical part of h following the same arguments used in particle transport. 
What is left we manipulate using tensor notation. We define a symmetric tensor 
R ( the "energy- weighted stress tensor", ( ]HH76P ). 



R= I d\ (mt;V2)vv/i^ 



The argument is finished exactly as before. The field terms in Eqn (A.67) will 
produce transport. 
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First we calculate the Fq part using the methods from particle transport. 
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What is left is 
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For the hgk part in Eqn( A.70), the integration by parts is more complex than 
with the analogous term in the particle transport equation. This is due to the 
f ^ factor. The v ■ V5A term is removed as with the density transport case of 



Eqn (A.62). The result is 
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We will later see that the second term of Eqn (A.71) will cancel with part of the 



electric field in Eqn (A. 65). The collisional terms of Eqn (A. 67) work similarly 
to density transport. 



(A.72) 



AV 



This covers all of the terms from Eqn (A. 67). The next step is to return to 



Eqn (A. 65) and calculate the electric field term. The scalar potential part of the 
electric field is an order larger than the ^ term, but as we will see most of it will 
cancel. This will result in the part of the scalar potential that does not average 
away acting on the same order as 
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(A.73) 



where we have integrated by parts in time and space between hnes one and two, 
between hnes two and three we use ^ = C(/, /) and, finally, we use gauss's 
law in velocity space and that all collisions conserve particles to obtain the final 
line. The second to last term that we obtain can be time averaged to zero. We 



find that the very last term will cancel the last term of Eqn (A. 71) when we 
plug in h for f in the equation. The Fq part of the final term is removed by 
integration and patch averaging to higher order. The remaining part of this final 
term is manipulated as follows (note that /i„c does not enter due to the small 
patch average of ip): 
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(A. 74) 
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This indeed cancels with the second term in Eqn (A. 71). Between hnes 1 and 2 



we used Eqn (A. 17). Between hnes 4 and 5 we integrated by parts in space to 



move the gradient onto the (p. AU that is left for us to calculate from the scalar 
potential term in E ■ v is 
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(A.75) 



We now substitute our solution for 6fi. The adiabatic/Boltzmann term gives 
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since the A\\ term is odd in v\\ and the A± gyro-averages to zero (or equivalently 
is odd in v±). This term will time average away, leading to the result that at 
this order the E ■ v^Fq term produces no heating. We collect what is left of the 



heating terms from Sfi of Eqn (A.75): 
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The final term we need to calculate is the collisional energy exchange terms 
between species are now included - note like particle collisions do not produce a 
loss of energy and thus do not appear (this is easy to prove). The interspecies 
collisional energy exchange is standard since to this order it is between Maxwellian 
species. Specifically we have 



d\ ( ^mV^Csrifs, fr) 
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where z/l*" is found in ( ]HS02p . It is ^fni^mi smaller than the ion collision rate 
which is itself ^/mjrni smaller than the electron collision rate. We now collect 



all the results given by Eqns (A.69), (A.71), (A.72), (A.77), and (A.78). 
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A. 8.1 Heat transport manipulation 



The purpose of this subsection is to manipulate Eqn (A.79) into a form that will 



readily apply to Sec. A. 6 We can rewrite the bracketed terms of Eqn (A.79) by 



using the gyro-kinetic equations (2.54) and (2.55). A change of variables trick 



will convert our integral into a gyro average with negligible error so that we may 
write 



AV 



+ 



T 
T 



gk 



2 dt 
1 
2 



Ay 



hgk {C{hgk))-^ 



AV 

(A.80) 



Ay 



2.55 



by ^ to obtain the right hand 



We have multiplied the gyrokinetic equation 
side of the equation. Most of this is zero. The first term time averages. The 
second term averages to zero over the annulus. (One uses that bo • Vf y = v_|_ ■ 
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(bo • Vbo) then gyro-averages.) Then the annulus average is done employing 
magnetic coordinates.) The third term annulus averages to higher order also 
using V ■ = and applying the divergence theorem. The fifth term is argued 
away by order by writing it as a full divergence with negligible error (V ■ is 
non-zero but small). Two terms remain: 



I d\ (q^h,,"^^^ = I d\ [v^ • VFo - {C{h,,,F^))]J (A.81) 

We use axi-symmetry and write Bo = V'i/'X V0+-F (■?/;) V0. (Note: F{ip) = i?0-Bo 
is not to be confused with Fq.) Then we integrate over the annulus to obtain 



where we have removed the Boltzman part of Fq by averaging over and inte- 
grating by parts in space to obtain a bo ■ V integral which annulus averages to 



higher order. Next we use the neoclassical equation (Eqn. |2.54) to write 
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= J rf^v^^(^ + .;|,bo-V/.„, + VD-VFo-(C(M)) 

= j d\(^ (VD ■ VFo - (C(/i„J))^ (A.83) 

where we have eliminated the first two terms on the second line by time average 
and annulus average respectively. Note that we only need to keep the F^ part 
of Fq here since patch averaging/gyro-averaging cleans up the rest. We need to 
rewrite the first term on the final line to show it is mostly collisional. We again 
use axi-symmetry to obtain 



186 



Ay 



mvlTR din Frr 
q dip 



■ bobo ■ 



(A.84) 



where we have integrated by parts between hne two and three (refer to the earher 
arguments in the particle transport derivation with pressure tensor P). Next, we 



again refer to the neoclassical equation (2.54) to write 



W||bo ■ Vhnc = {C{hnc)) - vd ■ VFo - g^'^ll^o ■ (A.85) 

The middle vd term will integrate to zero in f y upon substitution into equation 
IA.84I We will soon see that the last term will cancel. Let us collect the recent 
results: 
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Finally we need to rewrite the last term on the left hand side of equation A. 25 
The result is 
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We substitute these last three results into our heat transport equation (A. 25) to 
obtain 
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